ปฏิกิริยาสัมพันธ์ระหว่างการรักษากับเวลา: สัมประสิทธิ์ตัวไหนคือผลของการรักษา
On this page
Read the English version
บทคัดย่อ
ในการทดลองแบบสุ่มที่วัดผลซ้ำหลายครั้ง แบบจำลองถดถอยที่มีตัวแปรกลุ่มการรักษา (arm) ครั้งที่วัด (visit) และปฏิกิริยาสัมพันธ์ (interaction) ของทั้งสองจะพิมพ์สัมประสิทธิ์ออกมาหลายตัว โดยไม่มีตัวใดถูกระบุว่าเป็นผลของการรักษา บทความนี้จับคู่สัมประสิทธิ์แต่ละตัวกับค่าเฉลี่ยรายกลุ่มรายครั้งที่วัดที่มันเปรียบเทียบ เมื่อมีสองครั้งที่วัด สัมประสิทธิ์ของกลุ่มคือความต่างระหว่างกลุ่มที่จุดเริ่มต้น (baseline) สัมประสิทธิ์ของครั้งที่วัดคือการเปลี่ยนแปลงในกลุ่มควบคุม และปฏิกิริยาสัมพันธ์คือผลต่างของผลต่าง (difference in differences) ได้แก่การเปลี่ยนแปลงในกลุ่มรักษาลบการเปลี่ยนแปลงในกลุ่มควบคุม เมื่อมีสี่ครั้งที่วัดและลงรหัสแต่ละครั้งเป็นหมวดของตัวเอง ซึ่งเรียกว่าเวลาแบบจัดกลุ่ม (categorical time) ความต่างระหว่างกลุ่มในครั้งหลังเท่ากับสัมประสิทธิ์ของกลุ่มบวกสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ของครั้งนั้น การถือว่าเวลาเป็นตัวเลขให้ผลต่างของความชันหนึ่งค่า และสมมติว่าเป็นเส้นตรง การทดลองจำลองในผู้ใหญ่ 300 คนแสดงการคำนวณ เปรียบเทียบวิธีใช้คะแนนที่จุดเริ่มต้นสามวิธี และสร้างแบบจำลองสหสัมพันธ์ระหว่างครั้งที่วัด บทความนี้สรุปว่าต้องระบุครั้งที่วัดควบคู่กับผลของการรักษา และอ่านผลเป็นผลรวมของสัมประสิทธิ์ และความต่างนี้ในแบบที่ปรับด้วยคะแนนที่จุดเริ่มต้นมักแม่นยำกว่า
ตารางถดถอยในที่ประชุมเขียนรายงานการทดลอง
คณะผู้เขียนกำลังร่างรายงานหลักของการทดลองแบบสุ่ม การทดลองนี้รับผู้ใหญ่ 300 คนที่มีคะแนนอาการเรื้อรัง โดยคะแนนต่ำหมายถึงอาการน้อยลง วัดคะแนนที่สัปดาห์ที่ 0, 4, 8 และ 12 โดยสัปดาห์ที่ 0 คือก่อนได้รับยาครั้งแรก การรักษานี้เป็นของสมมติและข้อมูลเป็นข้อมูลจำลอง
ตารางถดถอยมีค่าจุดตัดแกน (intercept) หนึ่งค่า สัมประสิทธิ์ของกลุ่มหนึ่งตัว สัมประสิทธิ์ของครั้งที่วัดสามตัว และสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ระหว่างกลุ่มกับครั้งที่วัดสามตัว
ปฏิกิริยาสัมพันธ์ระหว่างการรักษากับเวลา (treatment-by-time interaction) ทำให้ความต่างระหว่างกลุ่มเปลี่ยนไปได้ในแต่ละครั้งที่วัด การทดสอบร่วม (joint test) ของสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ทั้งสามตัว ซึ่งเป็นการทดสอบครั้งเดียวว่าความต่างระหว่างกลุ่มเท่ากันทุกครั้งที่วัดหรือไม่ มีนัยสำคัญชัดเจน แพทย์ท่านหนึ่งถามว่า ในตัวเลขแปดตัวนี้ ตัวไหนคือผลของการรักษา
คำตอบประโยคเดียวคือ ผลของการรักษาที่สัปดาห์ที่ 12 เท่ากับสัมประสิทธิ์ของกลุ่มบวกสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ของสัปดาห์ที่ 12 และรายงานควรระบุสัปดาห์ที่ 12 เป็นครั้งที่วัดที่สนใจ
เนื่องจากการทดลองวัดคะแนนที่จุดเริ่มต้นก่อนการรักษา ความต่างนี้ในแบบที่ปรับด้วยคะแนนที่จุดเริ่มต้น (จุดเริ่มต้นร่วมแบบบังคับ หรือ ANCOVA ซึ่งนิยามไว้ในหัวข้อ 'สามวิธีใช้คะแนนสัปดาห์ที่ 0') มักเป็นค่าประมาณหลักที่แม่นยำกว่า ส่วนที่เหลือของบทความแสดงเหตุผล
สองกลุ่ม สองครั้งที่วัด: ค่าเฉลี่ยสี่ค่าและสัมประสิทธิ์สี่ตัว
เริ่มจากครั้งที่วัดเริ่มต้นและครั้งติดตามผลหนึ่งครั้ง ให้ $Y_{ij}$ เป็นคะแนนของผู้เข้าร่วมคนที่ $i$ ในครั้งที่วัด $j$ ให้ $\mathrm{arm}_i$ เป็น 1 สำหรับกลุ่มรักษาและ 0 สำหรับกลุ่มควบคุม และ $\mathrm{post}_j$ เป็น 1 ที่ครั้งติดตามผลและ 0 ที่จุดเริ่มต้น แบบจำลองของค่าเฉลี่ยคือ
$$\mathrm{E}(Y_{ij}) = \beta_0 + \beta_1\,\mathrm{arm}_i + \beta_2\,\mathrm{post}_j + \beta_3\,\mathrm{arm}_i\,\mathrm{post}_j$$ในสมการนี้ $\mathrm{E}(Y_{ij})$ คือคะแนนเฉลี่ยของกลุ่มและครั้งที่วัดที่กำหนด และ $\beta$ แต่ละตัวคือสัมประสิทธิ์ที่ต้องประมาณ เมื่อกำหนดตัวบ่งชี้ทั้งสองเป็น 0 หรือ 1 จะได้ค่าเฉลี่ยของเซลล์สี่ค่าในตาราง
ค่าเฉลี่ยของเซลล์สี่ค่าในแบบจำลองสองครั้งที่วัด
| กลุ่ม | จุดเริ่มต้น ($\mathrm{post}_j = 0$) | ติดตามผล ($\mathrm{post}_j = 1$) | การเปลี่ยนแปลง ติดตามผลลบจุดเริ่มต้น |
|---|---|---|---|
| กลุ่มควบคุม ($\mathrm{arm}_i = 0$) | $\beta_0$ | $\beta_0 + \beta_2$ | $\beta_2$ |
| กลุ่มรักษา ($\mathrm{arm}_i = 1$) | $\beta_0 + \beta_1$ | $\beta_0 + \beta_1 + \beta_2 + \beta_3$ | $\beta_2 + \beta_3$ |
| กลุ่มรักษาลบกลุ่มควบคุม | $\beta_1$ | $\beta_1 + \beta_3$ | $\beta_3$ |
อ่านเซลล์ทั้งสี่
สัมประสิทธิ์แต่ละตัวตอนนี้มีความหมายชัดเจน $\beta_0$ คือค่าเฉลี่ยของกลุ่มควบคุมที่จุดเริ่มต้น (baseline) ซึ่งเป็นครั้งที่วัดก่อนการรักษา และ $\beta_1$ คือความต่างระหว่างกลุ่มที่จุดเริ่มต้น $\beta_2$ คือการเปลี่ยนแปลงในกลุ่มควบคุม $\beta_3$ คือการเปลี่ยนแปลงในกลุ่มรักษาลบการเปลี่ยนแปลงในกลุ่มควบคุม หรือ ผลต่างของผลต่าง (difference in differences)
ความต่างระหว่างกลุ่มที่ครั้งติดตามผลคือ $\beta_1 + \beta_3$ ซึ่งเป็นผลรวม ไม่ใช่สัมประสิทธิ์ตัวใดตัวหนึ่งตามลำพัง คอนทราสต์ (contrast) คือผลรวมหรือผลต่างของสัมประสิทธิ์ที่ตอบคำถามเดียว เช่นความต่างระหว่างกลุ่มที่สัปดาห์ที่ 12
ระดับอ้างอิง (reference level) ของตัวแปรคือหมวดที่หมวดอื่นถูกนำมาเปรียบเทียบด้วย ในที่นี้คือหมวดที่ลงรหัสเป็น 0 ในแบบจำลองที่มีปฏิกิริยาสัมพันธ์ สัมประสิทธิ์ของกลุ่มหรือของครั้งที่วัดแต่ละตัวอ่านโดยให้ตัวแปรอีกตัวอยู่ที่ระดับอ้างอิงของมัน ในที่นี้คือจุดเริ่มต้นและกลุ่มควบคุม
ในการทดลอง จุดเริ่มต้นเกิดก่อนการรักษา การสุ่มจึงทำให้ความต่างที่คาดไว้ที่จุดเริ่มต้นเป็นศูนย์ ค่าประมาณของ $\beta_1$ จึงสะท้อนความไม่สมดุลที่เกิดจากความบังเอิญ ไม่ใช่ผลของการรักษา หากกำหนดให้ครั้งติดตามผลเป็นระดับอ้างอิง $\beta_1$ จะกลายเป็นความต่างที่ครั้งติดตามผล โดยค่าเฉลี่ยที่ประมาณได้ไม่เปลี่ยน
ตัวอย่างคำนวณด้วยมือ: ค่าเฉลี่ยจริงสี่ค่า สัมประสิทธิ์สี่ตัว
ตัวอย่างคำนวณด้วยมือ ในการทดลองจำลองที่ใช้ตลอดบทความนี้ ค่าเฉลี่ยจริงของคะแนน (ค่าจริงในข้อมูลจำลอง) ที่สัปดาห์ที่ 0 เท่ากับ 50 ในทั้งสองกลุ่ม ที่สัปดาห์ที่ 12 ค่าเฉลี่ยจริงเท่ากับ 44 ในกลุ่มควบคุม และ 38 ในกลุ่มรักษา ให้สัปดาห์ที่ 0 เป็นจุดเริ่มต้นและสัปดาห์ที่ 12 เป็นครั้งติดตามผล
-
ค่าเฉลี่ยกลุ่มควบคุมที่จุดเริ่มต้น
\[ \beta_0 = 50 \]
ค่าจุดตัดแกนคือกลุ่มควบคุมที่สัปดาห์ที่ 0
-
ความต่างที่จุดเริ่มต้น
\[ \beta_1 = 50 - 50 = 0 \]
การสุ่มทำให้ความต่างจริงที่จุดเริ่มต้นเป็นศูนย์
-
การเปลี่ยนแปลงในกลุ่มควบคุม
\[ \beta_2 = 44 - 50 = -6 \]
กลุ่มควบคุมดีขึ้น 6 คะแนนโดยไม่ได้รับการรักษา
-
ผลต่างของผลต่าง
\[ \beta_3 = (38 - 50) - (44 - 50) = -12 - (-6) = -6 \]
กลุ่มรักษาเปลี่ยนไป -12 คะแนน เทียบกับ -6 คะแนนในกลุ่มควบคุม
-
ผลของการรักษาที่สัปดาห์ที่ 12
\[ \beta_1 + \beta_3 = 0 + (-6) = -6 \]
ได้ตัวเลขเดียวกันโดยตรงจากค่าเฉลี่ยของสัปดาห์ที่ 12 คือ 38 ลบ 44 เท่ากับ -6
ผลลัพธ์: ความต่างจริงของคะแนนเฉลี่ยระหว่างกลุ่มที่สัปดาห์ที่ 12 คือ -6 คะแนน ซึ่งเท่ากับ $\beta_3$ ในที่นี้เพียงเพราะ $\beta_1 = 0$
สี่ครั้งที่วัด: สัมประสิทธิ์ปฏิกิริยาสัมพันธ์หนึ่งตัวต่อหนึ่งครั้ง
เมื่อวัดที่สัปดาห์ที่ 0, 4, 8 และ 12 แนวทางหนึ่งคือ เวลาแบบจัดกลุ่ม (categorical time) คือมีตัวบ่งชี้หนึ่งตัวต่อหนึ่งครั้งที่วัดหลังจุดเริ่มต้น โดยไม่สมมติรูปร่างใด ๆ ให้ $d_{kj}$ เป็น 1 เมื่อครั้งที่วัด $j$ เป็นครั้งที่ $k$ หลังจุดเริ่มต้น และเป็น 0 ในกรณีอื่น ดังนั้น $k$ = 1, 2 และ 3 หมายถึงสัปดาห์ที่ 4, 8 และ 12 แบบจำลองของค่าเฉลี่ยจึงเป็น
$$\mathrm{E}(Y_{ij}) = \beta_0 + \beta_1\,\mathrm{arm}_i + \sum_{k=1}^{3} \beta_{2,k}\,d_{kj} + \sum_{k=1}^{3} \beta_{3,k}\,\mathrm{arm}_i\,d_{kj}$$ครั้งที่วัดหลังแต่ละครั้งมีการเปลี่ยนแปลงในกลุ่มควบคุมของตัวเอง คือ $\beta_{2,k}$ และมีสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ของตัวเอง คือ $\beta_{3,k}$ ความต่างระหว่างกลุ่มที่ครั้งที่ $k$ คือ $\beta_1 + \beta_{3,k}$ เป็นผลรวมแบบเดิม เพียงแต่ดูทีละครั้ง เมื่อมีสัมประสิทธิ์แปดตัวสำหรับค่าเฉลี่ยของเซลล์แปดค่า แบบจำลองนี้อิ่มตัว (saturated) และถ้าไม่มีครั้งที่วัดใดขาดหาย มันจะทำซ้ำค่าเฉลี่ยของเซลล์ที่สังเกตได้อย่างตรงทุกค่า
การทดสอบร่วมของสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ทั้งสามถามว่า ความต่างระหว่างกลุ่มเท่ากันทุกครั้งที่วัดหรือไม่ แต่ผู้อ่านมักต้องการคำตอบอีกแบบ คือความต่างที่สัปดาห์ที่ 12 มีขนาดเท่าไร
ใช้เส้นตรงแทน: ผลต่างของความชันหนึ่งค่า
แนวทางที่สองถือว่าเวลาเป็นตัวเลข เรียกว่าเวลาเชิงเส้น (linear time) ให้ $t_j$ เป็นสัปดาห์ของครั้งที่วัด $j$ แบบจำลองเวลาเชิงเส้นคือ
$$\mathrm{E}(Y_{ij}) = \gamma_0 + \gamma_1\,\mathrm{arm}_i + \gamma_2\,t_j + \gamma_3\,\mathrm{arm}_i\,t_j$$ในสมการนี้ $\gamma_0$ และ $\gamma_1$ ทำหน้าที่เหมือน $\beta_0$ และ $\beta_1$ ส่วน $\gamma_2$ คือความชันของกลุ่มควบคุมในหน่วยคะแนนต่อสัปดาห์ และ $\gamma_3$ คือผลต่างของความชันระหว่างกลุ่ม ความต่างระหว่างกลุ่มที่สัปดาห์ $t$ คือ $\gamma_1 + \gamma_3 t$ สัมประสิทธิ์ตัวเดียวจึงบรรจุเส้นทางตามเวลาทั้งหมด ราคาที่ต้องจ่ายคือข้อสมมติว่าแต่ละกลุ่มเดินตามเส้นตรง
ในการทดลองจำลอง ความต่างจริงระหว่างกลุ่มคือ -2, -4 และ -6 คะแนนที่สัปดาห์ที่ 4, 8 และ 12 ซึ่งอยู่บนเส้นตรงเดียวกัน การลงรหัสทั้งสองแบบจึงมุ่งไปยังค่าเดียวกัน เมื่อแนวโน้มโค้ง เช่นตอบสนองเร็วแล้วราบลง แบบจำลองเชิงเส้นอาจระบุความต่างที่ครั้งสุดท้ายผิดไป เวลาแบบจัดกลุ่มไม่สมมติรูปร่างใด และเหมาะกับการทดลองที่มีครั้งที่วัดคงที่ไม่กี่ครั้ง [1]
เวลาแบบจัดกลุ่มเทียบกับเวลาเชิงเส้น
| ลักษณะ | เวลาแบบจัดกลุ่ม | เวลาเชิงเส้น |
|---|---|---|
| ตัวแปรเวลา | ตัวบ่งชี้หนึ่งตัวต่อหนึ่งครั้งที่วัดหลังจุดเริ่มต้น | สัปดาห์เป็นตัวเลข |
| พจน์ปฏิกิริยาสัมพันธ์ | หนึ่งตัวต่อหนึ่งครั้งที่วัดหลัง คือ $\beta_{3,k}$ | ผลต่างของความชันหนึ่งค่า คือ $\gamma_3$ |
| ความต่างระหว่างกลุ่มที่ครั้งที่วัดหนึ่ง | $\beta_1 + \beta_{3,k}$ | $\gamma_1 + \gamma_3 t$ |
| รูปร่างที่สมมติ | ไม่มี | เส้นตรงในแต่ละกลุ่ม |
| ความเสี่ยงหลัก | สัมประสิทธิ์มากขึ้น ช่วงเชื่อมั่นกว้างขึ้นเมื่อมีครั้งที่วัดมาก | เกิดความเอนเอียงที่ครั้งสุดท้ายเมื่อแนวโน้มโค้ง |
การทดลองจำลอง: แบบจำลองผสมให้ผลอย่างไร
ส่วนที่เหลือของบทความวิเคราะห์ชุดข้อมูลจำลองชุดเดียว ชุดนี้มีผู้เข้าร่วม 300 คน กลุ่มละ 150 คน แต่ละคนถูกวัดคะแนนสี่ครั้ง รวมเป็น 1,200 แถว โดยหนึ่งแถวคือผู้เข้าร่วมหนึ่งคนในหนึ่งครั้งที่วัด ค่าเฉลี่ยจริงตรงกับตัวอย่างคำนวณด้วยมือ ความต่างจริงระหว่างกลุ่มคือ -2, -4 และ -6 คะแนนที่สัปดาห์ที่ 4, 8 และ 12
คะแนนเฉลี่ยที่สังเกตได้ แยกตามกลุ่มและสัปดาห์
| กลุ่ม | สัปดาห์ที่ 0 (คะแนน) | สัปดาห์ที่ 4 (คะแนน) | สัปดาห์ที่ 8 (คะแนน) | สัปดาห์ที่ 12 (คะแนน) |
|---|---|---|---|---|
| กลุ่มควบคุม | 50.16 | 48.18 | 47.11 | 44.96 |
| กลุ่มรักษา | 49.92 | 46.18 | 41.88 | 37.54 |
ประมาณเวลาแบบจัดกลุ่มพร้อมผลสุ่ม (random effects)
ข้อมูลรายผู้เข้าร่วมไม่ได้เผยแพร่ ตารางค่าเฉลี่ยด้านบนเพียงพอสำหรับคำนวณสัมประสิทธิ์ด้วยมือ และโค้ดแสดงวิธีประมาณแบบจำลอง
แบบจำลองที่ใช้คือแบบจำลองผสมเชิงเส้น (linear mixed model) ซึ่งประมาณด้วย REML (ภาวะน่าจะเป็นสูงสุดแบบจำกัด หรือ restricted maximum likelihood ซึ่งเป็นวิธีปกติในการประมาณแบบจำลองเหล่านี้) ได้แก่แบบจำลองค่าเฉลี่ยของเวลาแบบจัดกลุ่มบวกพจน์ระดับผู้เข้าร่วมสองพจน์ ซึ่งเรียกว่าผลสุ่ม (random effects) จุดตัดแกนสุ่ม (random intercept) เลื่อนเส้นทางทั้งเส้นของผู้เข้าร่วมคนหนึ่งขึ้นหรือลง ความชันสุ่ม (random slope) เอียงเส้นนั้นในหน่วยคะแนนต่อสัปดาห์ พจน์เหล่านี้จำลองว่าคะแนนของผู้เข้าร่วมคนเดียวกันสัมพันธ์กันอย่างไร และไม่เปลี่ยนความหมายของสัมประสิทธิ์ค่าเฉลี่ย
เครื่องหมายหมวก (hat) แสดงค่าประมาณจากตัวอย่าง $\hat\beta_1$ เท่ากับ -0.24 คะแนน (ช่วงเชื่อมั่น 95% คือ -1.90 ถึง 1.42) ซึ่งเป็นความต่างที่จุดเริ่มต้นอันเกิดจากความบังเอิญ ปฏิกิริยาสัมพันธ์ของสัปดาห์ที่ 12 คือ $\hat\beta_{3,3}$ เท่ากับ -7.18 คะแนน (ช่วงเชื่อมั่น 95% คือ -8.98 ถึง -5.38)
ความต่างระหว่างกลุ่มที่สัปดาห์ที่ 12 คือ $\hat\beta_1 + \hat\beta_{3,3}$ เท่ากับ -7.42 คะแนน ตัวเลขสุดท้ายนี้เท่ากับ 37.54 ลบ 44.96 จากตาราง เพราะแบบจำลองค่าเฉลี่ยอิ่มตัว
โค้ดด้านล่างยังประมาณวิธีใช้คะแนนสัปดาห์ที่ 0 สามวิธี (ดูหัวข้อ 'สามวิธีใช้คะแนนสัปดาห์ที่ 0') และความแปรปรวนร่วมแบบไม่มีโครงสร้าง (unstructured covariance) ระหว่างจุดตัดแกนสุ่มกับความชันสุ่ม (ดูหัวข้อ 'ความสัมพันธ์ระหว่างครั้งที่วัด และสิ่งที่เปลี่ยนไปเมื่อมีครั้งที่วัดขาดหาย')
ความต่างระหว่างกลุ่มในแต่ละครั้งที่วัด: ค่าจริงและค่าประมาณ
| สัปดาห์ | ความต่างจริง (คะแนน) | ความต่างที่ประมาณได้ (คะแนน) | ช่วงเชื่อมั่น 95% (คะแนน) |
|---|---|---|---|
| 0 | 0 | -0.24 | -1.90 ถึง 1.42 |
| 4 | -2 | -2.00 | -3.67 ถึง -0.33 |
| 8 | -4 | -5.23 | -7.02 ถึง -3.44 |
| 12 | -6 | -7.42 | -9.41 ถึง -5.43 |
Stata: แบบจำลองผสม คอนทราสต์รายครั้งที่วัด และวิธีใช้คะแนนที่จุดเริ่มต้นสำหรับสัปดาห์ที่ 12
* Simulated trial: 300 participants, 1:1, symptom score at weeks 0, 4, 8, 12 (lower is better).
* Simulated data only: nothing here is evidence about any real drug or patient.
version 18
clear all
set more off
set linesize 120
* Load the simulated data (long format: one row per participant-visit).
import delimited using trial.csv, clear asdouble
label define armlab 0 "control" 1 "treatment"
label values arm armlab
* ---------------------------------------------------------------- mixed model, arm x visit
* Categorical visit, arm-by-visit interaction, random intercept and random slope on week (REML).
* stddeviations prints the random effects as standard deviations and their correlation.
mixed score i.arm##i.visit || id: week, covariance(unstructured) reml stddeviations
* beta1 = 1.arm: the difference between arms at week 0, the reference visit.
lincom 1.arm
* beta3,3 = 1.arm#3.visit: the week-12 interaction, a difference in differences.
lincom 1.arm#3.visit
* Difference between arms at each visit, beta1 + beta3,k (k = 1, 2, 3 for weeks 4, 8, 12).
contrast r.arm@visit, effects
forvalues k = 1/3 {
lincom 1.arm + 1.arm#`k'.visit
}
* Joint Wald test that the three interaction coefficients are zero.
testparm 1.arm#i.visit
* ---------------------------------------------------------------- constrained baseline
* Same random-effects structure, but one common baseline mean for both arms.
generate byte trt_wk4 = arm * (visit == 1)
generate byte trt_wk8 = arm * (visit == 2)
generate byte trt_wk12 = arm * (visit == 3)
mixed score i.visit trt_wk4 trt_wk8 trt_wk12 || id: week, covariance(unstructured) reml stddeviations
lincom trt_wk12
* ---------------------------------------------------------------- week-12 ANCOVA and change score
* Baseline score copied to every row of the participant.
bysort id (visit): generate double score0 = score[1]
generate double change = score - score0
* ANCOVA: week-12 score on arm, adjusted for the baseline score (OLS).
regress score i.arm score0 if visit == 3
* Change score: (week-12 minus week-0) on arm, no baseline adjustment (OLS).
regress change i.arm if visit == 3. * Simulated trial: 300 participants, 1:1, symptom score at weeks 0, 4, 8, 12
> (lower is better).
. * Simulated data only: nothing here is evidence about any real drug or patien
> t.
. version 18
. clear all
. set more off
. set linesize 120
.
. * Load the simulated data (long format: one row per participant-visit).
. import delimited using trial.csv, clear asdouble
(encoding automatically selected: ISO-8859-1)
(6 vars, 1,200 obs)
. label define armlab 0 "control" 1 "treatment"
. label values arm armlab
.
. * ---------------------------------------------------------------- mixed model, arm x visit
. * Categorical visit, arm-by-visit interaction, random intercept and random slope on week (REML).
. * stddeviations prints the random effects as standard deviations and their correlation.
. mixed score i.arm##i.visit || id: week, covariance(unstructured) reml stddeviations
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log restricted-likelihood = -3829.3213
Iteration 1: Log restricted-likelihood = -3829.3199
Iteration 2: Log restricted-likelihood = -3829.3199
Computing standard errors ...
Mixed-effects REML regression Number of obs = 1,200
Group variable: id Number of groups = 300
Obs per group:
min = 4
avg = 4.0
max = 4
Wald chi2(7) = 460.58
Log restricted-likelihood = -3829.3199 Prob > chi2 = 0.0000
----------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
arm |
treatment | -.2397333 .8462269 -0.28 0.777 -1.898308 1.418841
|
visit |
1 | -1.9842 .4868589 -4.08 0.000 -2.938426 -1.029974
2 | -3.056333 .5534852 -5.52 0.000 -4.141144 -1.971522
3 | -5.208 .649515 -8.02 0.000 -6.481026 -3.934974
|
arm#visit |
treatment#1 | -1.761667 .6885225 -2.56 0.011 -3.111146 -.4121874
treatment#2 | -4.987867 .7827463 -6.37 0.000 -6.522021 -3.453712
treatment#3 | -7.177467 .9185529 -7.81 0.000 -8.977797 -5.377136
|
_cons | 50.16413 .5983728 83.83 0.000 48.99134 51.33692
----------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(week) | .4654107 .0387788 .3952873 .5479739
sd(_cons) | 6.137018 .3306216 5.52205 6.820472
corr(week,_cons) | -.1110897 .092756 -.2872992 .0723929
-----------------------------+------------------------------------------------
sd(Residual) | 4.00556 .1160179 3.784503 4.23953
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 684.31 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
. * beta1 = 1.arm: the difference between arms at week 0, the reference visit.
. lincom 1.arm
( 1) [score]1.arm = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -.2397333 .8462269 -0.28 0.777 -1.898308 1.418841
------------------------------------------------------------------------------
. * beta3,3 = 1.arm#3.visit: the week-12 interaction, a difference in differences.
. lincom 1.arm#3.visit
( 1) [score]1.arm#3.visit = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -7.177467 .9185529 -7.81 0.000 -8.977797 -5.377136
------------------------------------------------------------------------------
. * Difference between arms at each visit, beta1 + beta3,k (k = 1, 2, 3 for weeks 4, 8, 12).
. contrast r.arm@visit, effects
Contrasts of marginal linear predictions
Margins: asbalanced
-------------------------------------------------------------
| df chi2 P>chi2
--------------------------+----------------------------------
score |
arm@visit |
(treatment vs control) 0 | 1 0.08 0.7769
(treatment vs control) 1 | 1 5.50 0.0190
(treatment vs control) 2 | 1 32.80 0.0000
(treatment vs control) 3 | 1 53.39 0.0000
Joint | 4 79.42 0.0000
-------------------------------------------------------------
-------------------------------------------------------------------------------------------
| Contrast Std. err. z P>|z| [95% conf. interval]
--------------------------+----------------------------------------------------------------
score |
arm@visit |
(treatment vs control) 0 | -.2397333 .8462269 -0.28 0.777 -1.898308 1.418841
(treatment vs control) 1 | -2.0014 .8535013 -2.34 0.019 -3.674232 -.3285682
(treatment vs control) 2 | -5.2276 .9128241 -5.73 0.000 -7.016702 -3.438498
(treatment vs control) 3 | -7.4172 1.015111 -7.31 0.000 -9.406781 -5.427619
-------------------------------------------------------------------------------------------
. forvalues k = 1/3 {
2. lincom 1.arm + 1.arm#`k'.visit
3. }
( 1) [score]1.arm + [score]1.arm#1.visit = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -2.0014 .8535013 -2.34 0.019 -3.674232 -.3285682
------------------------------------------------------------------------------
( 1) [score]1.arm + [score]1.arm#2.visit = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -5.2276 .9128241 -5.73 0.000 -7.016702 -3.438498
------------------------------------------------------------------------------
( 1) [score]1.arm + [score]1.arm#3.visit = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -7.4172 1.015111 -7.31 0.000 -9.406781 -5.427619
------------------------------------------------------------------------------
. * Joint Wald test that the three interaction coefficients are zero.
. testparm 1.arm#i.visit
( 1) [score]1.arm#1.visit = 0
( 2) [score]1.arm#2.visit = 0
( 3) [score]1.arm#3.visit = 0
chi2( 3) = 70.56
Prob > chi2 = 0.0000
.
. * ---------------------------------------------------------------- constrained baseline
. * Same random-effects structure, but one common baseline mean for both arms.
. generate byte trt_wk4 = arm * (visit == 1)
. generate byte trt_wk8 = arm * (visit == 2)
. generate byte trt_wk12 = arm * (visit == 3)
. mixed score i.visit trt_wk4 trt_wk8 trt_wk12 || id: week, covariance(unstructured) reml stddeviations
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log restricted-likelihood = -3830.1128
Iteration 1: Log restricted-likelihood = -3830.1114
Iteration 2: Log restricted-likelihood = -3830.1114
Computing standard errors ...
Mixed-effects REML regression Number of obs = 1,200
Group variable: id Number of groups = 300
Obs per group:
min = 4
avg = 4.0
max = 4
Wald chi2(6) = 460.62
Log restricted-likelihood = -3830.1114 Prob > chi2 = 0.0000
----------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
visit |
1 | -1.945501 .4673096 -4.16 0.000 -2.861411 -1.029591
2 | -3.014831 .5337274 -5.65 0.000 -4.060917 -1.968744
3 | -5.163694 .6303518 -8.19 0.000 -6.399161 -3.928227
|
trt_wk4 | -1.839065 .6320851 -2.91 0.004 -3.077929 -.6002009
trt_wk8 | -5.070871 .7258901 -6.99 0.000 -6.49359 -3.648153
trt_wk12 | -7.266078 .8636518 -8.41 0.000 -8.958804 -5.573352
_cons | 50.04427 .4225711 118.43 0.000 49.21604 50.87249
----------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
id: Unstructured |
sd(week) | .4652983 .0387653 .3951986 .5478321
sd(_cons) | 6.125979 .3297757 5.51256 6.807657
corr(week,_cons) | -.1098813 .0927455 -.2861117 .0735397
-----------------------------+------------------------------------------------
sd(Residual) | 4.005283 .1159944 3.78427 4.239204
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 684.45 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.
. lincom trt_wk12
( 1) [score]trt_wk12 = 0
------------------------------------------------------------------------------
score | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
(1) | -7.266078 .8636518 -8.41 0.000 -8.958804 -5.573352
------------------------------------------------------------------------------
.
. * ---------------------------------------------------------------- week-12 ANCOVA and change score
. * Baseline score copied to every row of the participant.
. bysort id (visit): generate double score0 = score[1]
. generate double change = score - score0
. * ANCOVA: week-12 score on arm, adjusted for the baseline score (OLS).
. regress score i.arm score0 if visit == 3
Source | SS df MS Number of obs = 300
-------------+---------------------------------- F(2, 297) = 96.75
Model | 11054.4924 2 5527.24619 Prob > F = 0.0000
Residual | 16967.1828 297 57.1285616 R-squared = 0.3945
-------------+---------------------------------- Adj R-squared = 0.3904
Total | 28021.6752 299 93.7179772 Root MSE = 7.5583
------------------------------------------------------------------------------
score | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
arm |
treatment | -7.25869 .8728811 -8.32 0.000 -8.976506 -5.540875
score0 | .6611917 .0600397 11.01 0.000 .5430346 .7793488
_cons | 11.78802 3.074414 3.83 0.000 5.737628 17.83842
------------------------------------------------------------------------------
. * Change score: (week-12 minus week-0) on arm, no baseline adjustment (OLS).
. regress change i.arm if visit == 3
Source | SS df MS Number of obs = 300
-------------+---------------------------------- F(1, 298) = 61.29
Model | 3863.70208 1 3863.70208 Prob > F = 0.0000
Residual | 18786.4013 298 63.0416152 R-squared = 0.1706
-------------+---------------------------------- Adj R-squared = 0.1678
Total | 22650.1034 299 75.7528542 Root MSE = 7.9399
------------------------------------------------------------------------------
change | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
arm |
treatment | -7.177467 .9168178 -7.83 0.000 -8.981724 -5.373209
_cons | -5.208 .6482881 -8.03 0.000 -6.483803 -3.932197
------------------------------------------------------------------------------R: แบบจำลองชุดเดียวกันด้วย lme4
# Simulated trial: 300 participants, 1:1, symptom score at weeks 0, 4, 8, 12 (lower is better).
# Simulated data only: nothing here is evidence about any real drug or patient.
suppressPackageStartupMessages(library(lme4))
z975 <- qnorm(0.975)
# Load the simulated data (long format: one row per participant-visit).
trial <- read.csv("trial.csv")
trial <- trial[order(trial$id, trial$visit), ]
trial$visitf <- factor(trial$visit)
# ---------------------------------------------------------------- mixed model, arm x visit
# Categorical visit, arm-by-visit interaction, random intercept and random slope on week (REML).
lmm <- lmer(score ~ arm * visitf + (week | id), data = trial, REML = TRUE)
print(summary(lmm), correlation = FALSE)
b <- fixef(lmm)
V <- as.matrix(vcov(lmm))
# Wald estimate, SE and 95% CI (normal) of a sum of coefficients.
lc <- function(terms) {
L <- setNames(numeric(length(b)), names(b))
L[terms] <- 1
est <- sum(L * b)
se <- sqrt(drop(t(L) %*% V %*% L))
c(estimate = est, se = se, lo = est - z975 * se, hi = est + z975 * se)
}
# beta1 = arm: the difference between arms at week 0, the reference visit.
print(round(lc("arm"), 4))
# beta3,3 = arm:visitf3: the week-12 interaction, a difference in differences.
print(round(lc("arm:visitf3"), 4))
# Difference between arms at each visit, beta1 + beta3,k (k = 1, 2, 3 for weeks 4, 8, 12).
diffs <- t(sapply(0:3, function(k) lc(if (k == 0) "arm" else c("arm", paste0("arm:visitf", k)))))
rownames(diffs) <- paste("week", 4 * (0:3))
print(round(diffs, 4))
# Joint Wald test that the three interaction coefficients are zero.
idx <- paste0("arm:visitf", 1:3)
chi2 <- drop(t(b[idx]) %*% solve(V[idx, idx]) %*% b[idx])
cat(sprintf("Joint test: chi2(3) = %.4f, p = %.4g\n", chi2, pchisq(chi2, 3, lower.tail = FALSE)))
# ---------------------------------------------------------------- constrained baseline
# Same random-effects structure, but one common baseline mean for both arms.
trial$trt_wk4 <- trial$arm * (trial$visit == 1)
trial$trt_wk8 <- trial$arm * (trial$visit == 2)
trial$trt_wk12 <- trial$arm * (trial$visit == 3)
cb <- lmer(score ~ visitf + trt_wk4 + trt_wk8 + trt_wk12 + (week | id), data = trial, REML = TRUE)
est <- fixef(cb)[["trt_wk12"]]
se <- sqrt(as.matrix(vcov(cb))["trt_wk12", "trt_wk12"])
print(round(c(estimate = est, se = se, lo = est - z975 * se, hi = est + z975 * se), 4))
# ---------------------------------------------------------------- week-12 ANCOVA and change score
# One row per participant: baseline and week-12 scores side by side.
wide <- data.frame(id = trial$id[trial$visit == 0], arm = trial$arm[trial$visit == 0],
score0 = trial$score[trial$visit == 0], score12 = trial$score[trial$visit == 3])
# ANCOVA: week-12 score on arm, adjusted for the baseline score (OLS).
anc <- lm(score12 ~ arm + score0, data = wide)
print(summary(anc))
print(round(confint(anc), 4))
# Change score: (week-12 minus week-0) on arm, no baseline adjustment (OLS).
chg <- lm(I(score12 - score0) ~ arm, data = wide)
print(summary(chg))
print(round(confint(chg), 4)) > suppressPackageStartupMessages(library(lme4))
> z975 <- qnorm(0.975)
> trial <- read.csv("trial.csv")
> trial <- trial[order(trial$id, trial$visit), ]
> trial$visitf <- factor(trial$visit)
> lmm <- lmer(score ~ arm * visitf + (week | id), data = trial,
+ REML = TRUE)
> print(summary(lmm), correlation = FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ arm * visitf + (week | id)
Data: trial
REML criterion at convergence: 7658.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.8884 -0.5185 -0.0007 0.5334 3.1753
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 37.6629 6.1370
week 0.2166 0.4654 -0.11
Residual 16.0445 4.0056
Number of obs: 1200, groups: id, 300
Fixed effects:
Estimate Std. Error t value
(Intercept) 50.1641 0.5984 83.834
arm -0.2397 0.8462 -0.283
visitf1 -1.9842 0.4869 -4.076
visitf2 -3.0563 0.5535 -5.522
visitf3 -5.2080 0.6495 -8.018
arm:visitf1 -1.7617 0.6885 -2.559
arm:visitf2 -4.9879 0.7827 -6.372
arm:visitf3 -7.1775 0.9186 -7.814
> b <- fixef(lmm)
> V <- as.matrix(vcov(lmm))
> lc <- function(terms) {
+ L <- setNames(numeric(length(b)), names(b))
+ L[terms] <- 1
+ est <- sum(L * b)
+ se <- sqrt(drop(t(L) %*% V %*% L))
+ c(estimate = est, se = se, lo = est - z975 * se, hi = est +
+ z975 * se)
+ }
> print(round(lc("arm"), 4))
estimate se lo hi
-0.2397 0.8462 -1.8983 1.4188
> print(round(lc("arm:visitf3"), 4))
estimate se lo hi
-7.1775 0.9186 -8.9778 -5.3771
> diffs <- t(sapply(0:3, function(k) lc(if (k == 0) "arm" else c("arm",
+ paste0("arm:visitf", k)))))
> rownames(diffs) <- paste("week", 4 * (0:3))
> print(round(diffs, 4))
estimate se lo hi
week 0 -0.2397 0.8462 -1.8983 1.4188
week 4 -2.0014 0.8535 -3.6742 -0.3286
week 8 -5.2276 0.9128 -7.0167 -3.4385
week 12 -7.4172 1.0151 -9.4068 -5.4276
> idx <- paste0("arm:visitf", 1:3)
> chi2 <- drop(t(b[idx]) %*% solve(V[idx, idx]) %*%
+ b[idx])
> cat(sprintf("Joint test: chi2(3) = %.4f, p = %.4g\n",
+ chi2, pchisq(chi2, 3, lower.tail = FALSE)))
Joint test: chi2(3) = 70.5554, p = 3.246e-15
> trial$trt_wk4 <- trial$arm * (trial$visit == 1)
> trial$trt_wk8 <- trial$arm * (trial$visit == 2)
> trial$trt_wk12 <- trial$arm * (trial$visit == 3)
> cb <- lmer(score ~ visitf + trt_wk4 + trt_wk8 + trt_wk12 +
+ (week | id), data = trial, REML = TRUE)
> est <- fixef(cb)[["trt_wk12"]]
> se <- sqrt(as.matrix(vcov(cb))["trt_wk12", "trt_wk12"])
> print(round(c(estimate = est, se = se, lo = est -
+ z975 * se, hi = est + z975 * se), 4))
estimate se lo hi
-7.2661 0.8636 -8.9588 -5.5734
> wide <- data.frame(id = trial$id[trial$visit == 0],
+ arm = trial$arm[trial$visit == 0], score0 = trial$score[trial$visit ==
+ 0], score12 = trial$score[trial$visit == 3])
> anc <- lm(score12 ~ arm + score0, data = wide)
> print(summary(anc))
Call:
lm(formula = score12 ~ arm + score0, data = wide)
Residuals:
Min 1Q Median 3Q Max
-24.6458 -5.1035 0.8545 4.7741 15.7149
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.78802 3.07441 3.834 0.000154 ***
arm -7.25869 0.87288 -8.316 3.32e-15 ***
score0 0.66119 0.06004 11.013 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.558 on 297 degrees of freedom
Multiple R-squared: 0.3945, Adjusted R-squared: 0.3904
F-statistic: 96.75 on 2 and 297 DF, p-value: < 2.2e-16
> print(round(confint(anc), 4))
2.5 % 97.5 %
(Intercept) 5.7376 17.8384
arm -8.9765 -5.5409
score0 0.5430 0.7793
> chg <- lm(I(score12 - score0) ~ arm, data = wide)
> print(summary(chg))
Call:
lm(formula = I(score12 - score0) ~ arm, data = wide)
Residuals:
Min 1Q Median 3Q Max
-25.0045 -5.0189 0.9755 5.0555 16.9255
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.2080 0.6483 -8.033 2.23e-14 ***
arm -7.1775 0.9168 -7.829 8.69e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.94 on 298 degrees of freedom
Multiple R-squared: 0.1706, Adjusted R-squared: 0.1678
F-statistic: 61.29 on 1 and 298 DF, p-value: 8.689e-14
> print(round(confint(chg), 4))
2.5 % 97.5 %
(Intercept) -6.4838 -3.9322
arm -8.9817 -5.3732 สามวิธีใช้คะแนนสัปดาห์ที่ 0
การทดลองที่วัดผลลัพธ์ที่จุดเริ่มต้นสามารถใช้คะแนนนั้นได้อย่างน้อยสามวิธี เมื่อเป้าหมายคือความต่างที่สัปดาห์ที่ 12 [2, 3]
- ANCOVA (analysis of covariance หรือการวิเคราะห์ความแปรปรวนร่วม): ถดถอยคะแนนสัปดาห์ที่ 12 บนกลุ่มและคะแนนที่จุดเริ่มต้น
- คะแนนการเปลี่ยนแปลง (change score): ถดถอยคะแนนสัปดาห์ที่ 12 ลบคะแนนที่จุดเริ่มต้น บนกลุ่ม
- จุดเริ่มต้นร่วมแบบบังคับ (constrained baseline): ประมาณแบบจำลองระยะยาวโดยให้ทั้งสองกลุ่มมีค่าเฉลี่ยที่จุดเริ่มต้นร่วมกันหนึ่งค่า เพราะการสุ่มเกิดหลังสัปดาห์ที่ 0
ในการทดลองแบบสุ่ม ทั้งสามวิธีมุ่งไปที่ความต่างที่สัปดาห์ที่ 12 ค่าเดียวกัน เพราะความต่างที่คาดไว้ที่จุดเริ่มต้นเป็นศูนย์ วิธีเหล่านี้ต่างกันที่การจัดการกับความต่างที่จุดเริ่มต้นซึ่งเกิดจากความบังเอิญ และจึงต่างกันที่ความแม่นยำ
ความต่างระหว่างกลุ่มที่สัปดาห์ที่ 12 แยกตามวิธีใช้คะแนนที่จุดเริ่มต้น
| วิธี | ค่าประมาณ (คะแนน) | ค่าคลาดเคลื่อนมาตรฐาน (คะแนน) | ช่วงเชื่อมั่น 95% (คะแนน) |
|---|---|---|---|
| แบบจำลองผสม คอนทราสต์ที่ไม่ปรับที่สัปดาห์ที่ 12 | -7.42 | 1.02 | -9.41 ถึง -5.43 |
| ANCOVA | -7.26 | 0.87 | -8.98 ถึง -5.54 |
| คะแนนการเปลี่ยนแปลง | -7.18 | 0.92 | -8.98 ถึง -5.37 |
| จุดเริ่มต้นร่วมแบบบังคับ | -7.27 | 0.86 | -8.96 ถึง -5.57 |
ค่าประมาณทั้งสามเชื่อมโยงกันอย่างไร
สำหรับข้อมูลครบถ้วนชุดนี้ ความเชื่อมโยงเป็นไปอย่างตรงทุกตัวเลข ค่าประมาณของคะแนนการเปลี่ยนแปลงคือ -7.18 เท่ากับ $\hat\beta_{3,3}$ เพราะทั้งสองคือการเปลี่ยนแปลงในกลุ่มรักษาลบการเปลี่ยนแปลงในกลุ่มควบคุม ช่วงเชื่อมั่นของคะแนนการเปลี่ยนแปลงต่างกันที่ทศนิยมตำแหน่งที่สอง เพราะได้จากการถดถอยแยกอีกชุดหนึ่ง
ค่าประมาณของ ANCOVA เท่ากับความต่างดิบที่สัปดาห์ที่ 12 ลบ 0.66 คูณความต่างที่จุดเริ่มต้น โดย 0.66 คือสัมประสิทธิ์ของคะแนนที่จุดเริ่มต้นที่ประมาณได้ ส่วนคะแนนการเปลี่ยนแปลงตรึงตัวคูณนั้นไว้ที่ 1
เพราะ ANCOVA ประมาณตัวคูณจากข้อมูล จึงมักแม่นยำที่สุดในบรรดาวิธีง่าย ๆ [2] ค่าประมาณจากวิธีจุดเริ่มต้นร่วมแบบบังคับคือ -7.27 ใกล้เคียงกับ ANCOVA และมีค่าคลาดเคลื่อนมาตรฐานใกล้กัน
คอนทราสต์ที่สัปดาห์ที่ 12 แบบไม่ปรับจากแบบจำลองผสมมีค่าคลาดเคลื่อนมาตรฐานมากที่สุดคือ 1.02 คะแนน และช่วงเชื่อมั่นกว้างที่สุด มันเปรียบเทียบค่าเฉลี่ยของสัปดาห์ที่ 12 โดยตรงและไม่ใช้ข้อมูลในคะแนนที่จุดเริ่มต้น ซึ่ง ANCOVA และจุดเริ่มต้นร่วมแบบบังคับใช้ตัดความแปรปรวนจากความบังเอิญออกบางส่วน
ความสัมพันธ์ระหว่างครั้งที่วัด และสิ่งที่เปลี่ยนไปเมื่อมีครั้งที่วัดขาดหาย
คะแนนของผู้เข้าร่วมคนเดียวกันสัมพันธ์กัน คนที่เริ่มสูงมักอยู่ในระดับสูงต่อไป หากถือว่า 1,200 แถวเป็นคนอิสระ 1,200 คน ค่าคลาดเคลื่อนมาตรฐานจะผิด แบบจำลองจึงต้องมีโครงสร้างความแปรปรวนร่วม (covariance structure) ของคะแนนข้ามครั้งที่วัด
แบบจำลองข้างต้นใช้จุดตัดแกนสุ่มและความชันสุ่มพร้อม ความแปรปรวนร่วมแบบไม่มีโครงสร้าง (unstructured covariance) ระหว่างทั้งสอง ความแปรปรวนและสหสัมพันธ์ของมันจึงถูกประมาณอย่างอิสระ ส่วนเบี่ยงเบนมาตรฐานที่ประมาณได้คือ 6.14 คะแนนสำหรับจุดตัดแกน และ 0.47 คะแนนต่อสัปดาห์สำหรับความชัน สหสัมพันธ์ระหว่างทั้งสองคือ -0.11 และส่วนเบี่ยงเบนมาตรฐานของส่วนเหลือ (residual) คือ 4.01 คะแนน
ทางเลือกที่ใช้กันแพร่หลายคือ MMRM หรือแบบจำลองผสมสำหรับการวัดซ้ำ (mixed model for repeated measures) ซึ่งคงเวลาแบบจัดกลุ่มไว้ ตัดผลสุ่มออก และใส่ความแปรปรวนร่วมแบบไม่มีโครงสร้างบนส่วนเหลือ แต่ละครั้งที่วัดมีความแปรปรวนของตัวเอง และแต่ละคู่ของครั้งที่วัดมีสหสัมพันธ์ของตัวเอง [1]
เมื่อข้อมูลครบถ้วนและแบบจำลองค่าเฉลี่ยอิ่มตัว ทั้งสองแบบให้ค่าประมาณจุดเท่ากัน และต่างกันเฉพาะค่าคลาดเคลื่อนมาตรฐาน เมื่อมีครั้งที่วัดขาดหาย ค่าประมาณก็เปลี่ยนด้วย แบบจำลองที่อิงภาวะน่าจะเป็น (likelihood-based model) เลือกค่าพารามิเตอร์ที่ทำให้คะแนนที่สังเกตได้มีโอกาสเกิดขึ้นสูงที่สุด แบบจำลองนี้ใช้ทุกครั้งที่วัดที่สังเกตได้ของผู้เข้าร่วมแต่ละคน จึงยืมข้อมูล (borrow strength) จากครั้งเหล่านั้น
การยืมข้อมูลจากครั้งที่วัดที่สังเกตได้ใช้ได้ภายใต้ข้อมูลขาดหายแบบสุ่ม (missing at random, MAR) คือการขาดหายของครั้งที่วัดขึ้นกับข้อมูลที่สังเกตได้แล้วเท่านั้น และยังต้องกำหนดแบบจำลองค่าเฉลี่ยและความแปรปรวนร่วมถูกต้องด้วย ดังนั้นจุดตัดแกนสุ่มและความชันสุ่มที่เข้ากับข้อมูลได้ไม่ดีอาจทำให้คอนทราสต์ที่สัปดาห์ที่ 12 เปลี่ยนไปเมื่อมีครั้งที่วัดขาดหาย นี่เป็นเหตุผลหนึ่งที่มักกำหนดความแปรปรวนร่วมแบบไม่มีโครงสร้างข้ามครั้งที่วัดไว้ล่วงหน้าเมื่อคาดว่าจะมีผู้ถอนตัว [1]
ระบุครั้งที่วัดไว้ใน estimand
Estimand (ปริมาณเป้าหมายของการประมาณ) คือปริมาณที่แน่ชัดซึ่งการศึกษาตั้งใจประมาณ ภาคผนวก ICH E9(R1) อธิบาย estimand ด้วยองค์ประกอบห้าประการ ได้แก่ การรักษา ประชากร ตัวแปร (ผลลัพธ์ที่วัด) การจัดการกับเหตุการณ์ระหว่างการศึกษา (intercurrent events) และมาตรวัดสรุประดับประชากร (population-level summary) [4] เหตุการณ์ระหว่างการศึกษาคือเหตุการณ์ที่เกิดหลังเริ่มการรักษา เช่นการหยุดการรักษาหรือการใช้ยาเพิ่มเติม ที่มีผลต่อการวัดผลลัพธ์ได้หรือไม่ หรือต่อวิธีอ่านผลลัพธ์
ในการทดลองที่วัดผลสี่ครั้ง "ผลของการรักษา" ยังไม่ใช่ estimand เพราะความต่างหลังจุดเริ่มต้นทั้งสามไม่จำเป็นต้องเท่ากัน ข้อความที่ใช้ได้คือ เปรียบเทียบการรักษากับกลุ่มควบคุมในประชากรที่ถูกสุ่ม ใช้คะแนนอาการที่สัปดาห์ที่ 12 เป็นตัวแปร และสรุปเป็นความต่างของคะแนนเฉลี่ยระหว่างกลุ่ม estimand ฉบับสมบูรณ์ยังระบุด้วยว่าจะจัดการกับเหตุการณ์ระหว่างการศึกษาอย่างไร
เมื่อระบุครั้งที่วัดแล้ว ค่าประมาณคือคอนทราสต์ที่ครั้งนั้น ได้แก่ $\hat\beta_1 + \hat\beta_{3,k}$ หรือสัมประสิทธิ์ของกลุ่มจาก ANCOVA ค่าเฉลี่ยตลอดสัปดาห์ที่ 4, 8 และ 12 เป็น estimand อีกแบบ ซึ่งคำนวณเป็นการรวมเชิงเส้น (linear combination) หมายถึงผลรวมถ่วงน้ำหนักของสัมประสิทธิ์ (ในที่นี้คือค่าเฉลี่ยของความต่างหลังจุดเริ่มต้นทั้งสามค่า พร้อมช่วงเชื่อมั่นของมันเอง)
ความเข้าใจผิดที่พบบ่อยและวิธีแก้
-
"beta1 คือผลของการรักษา"
เมื่อมีปฏิกิริยาสัมพันธ์ สัมประสิทธิ์ของกลุ่มอ่านที่ครั้งที่วัดอ้างอิง ที่จุดเริ่มต้นก่อนการรักษา มันวัดความไม่สมดุลที่เกิดจากความบังเอิญ
วิธีแก้: ในแบบจำลองที่มีกลุ่มการรักษา ครั้งที่วัด และปฏิกิริยาสัมพันธ์ระหว่างทั้งสอง beta1 คือผลต่างระหว่างกลุ่ม ณ ครั้งที่วัดอ้างอิง ผลของการรักษา ณ ครั้งที่วัดหลังจากนั้นคือ beta1 บวกสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ของครั้งนั้น และ estimand ต้องระบุครั้งที่วัดเสมอ
-
"beta3 คือผลของการรักษาเสมอ"
ในการทดลองแบบสุ่มที่วัดคะแนนที่จุดเริ่มต้นก่อนการรักษา $\hat\beta_{3,k}$ คือค่าประมาณแบบคะแนนการเปลี่ยนแปลงของความต่างที่สัปดาห์ $k$ ค่าเดียวกัน คือ -7.18 คะแนนที่สัปดาห์ที่ 12 เทียบกับ -7.42 จาก $\hat\beta_1 + \hat\beta_{3,3}$ ทั้งสองค่าต่างกันเพียงเพราะความไม่สมดุลที่จุดเริ่มต้นซึ่งเกิดจากความบังเอิญ และ ANCOVA (-7.26) หรือจุดเริ่มต้นร่วมแบบบังคับ (-7.27) มักแม่นยำกว่า การอ่าน $\beta_{3,k}$ ตัวเดียวเป็นผลของการรักษาจะผิดเมื่อกลุ่มไม่ได้มาจากการสุ่ม เมื่อครั้งที่วัดอ้างอิงอยู่หลังเริ่มการรักษา หรือเมื่อไม่ได้วัดค่าเริ่มต้นก่อนการรักษา
วิธีแก้: เมื่อมีครั้งที่วัดเริ่มต้น ผลของการรักษาที่ครั้งหลังคือผลรวม $\beta_1 + \beta_{3,k}$ ซึ่งเท่ากับ $\beta_{3,k}$ ก็ต่อเมื่อ $\beta_1 = 0$
-
"สัมประสิทธิ์ของกลุ่มคือผลเฉลี่ยตลอดทุกครั้งที่วัด"
เมื่อมีปฏิกิริยาสัมพันธ์ในแบบจำลอง ไม่มีสัมประสิทธิ์ตัวเดียวที่ใช้ได้กับทุกครั้งที่วัด
วิธีแก้: ระบุค่าเฉลี่ยนั้นเป็น estimand ของมันเอง และคำนวณเป็นการรวมเชิงเส้นพร้อมช่วงเชื่อมั่นของมันเอง
-
"สัมประสิทธิ์ของกลุ่มคือผลที่ครั้งสุดท้าย"
จริงเฉพาะเมื่อครั้งสุดท้ายเป็นระดับอ้างอิง การเปลี่ยนรหัสเปลี่ยนสัมประสิทธิ์ แต่ไม่เปลี่ยนค่าเฉลี่ยที่ประมาณได้
วิธีแก้: ตรวจครั้งที่วัดอ้างอิง หรือรายงานคอนทราสต์ที่ครั้งที่วัดซึ่งระบุไว้ ซึ่งไม่ขึ้นกับการลงรหัส
-
"สัมประสิทธิ์ของครั้งที่วัดแสดงว่าการรักษาออกฤทธิ์อย่างไรตามเวลา"
$\beta_{2,k}$ แต่ละตัวคือการเปลี่ยนแปลงจากจุดเริ่มต้นในกลุ่มควบคุมเท่านั้น
วิธีแก้: การเปลี่ยนแปลงในกลุ่มรักษาคือ $\beta_{2,k} + \beta_{3,k}$
-
"ตัดปฏิกิริยาสัมพันธ์ออกแล้วจะได้ผลของการรักษาร่วมหนึ่งค่า"
เมื่อแบบจำลองมีครั้งที่วัดเริ่มต้นก่อนการรักษา การตัดนี้บังคับให้ความต่างที่สัปดาห์ที่ 0 เท่ากับความต่างทุกครั้งหลัง สัมประสิทธิ์ตัวเดียวจึงผสมความต่างที่จุดเริ่มต้น ซึ่งเป็นศูนย์โดยการออกแบบ เข้ากับความต่างครั้งหลัง โดยมักถูกดึงเข้าหาศูนย์
วิธีแก้: คงปฏิกิริยาสัมพันธ์ไว้ หรือปรับด้วยคะแนนที่จุดเริ่มต้นและสร้างแบบจำลองเฉพาะครั้งที่วัดหลังจุดเริ่มต้น โดยระบุผลร่วมใด ๆ เป็นข้อสมมติ
-
"แบบจำลองความแปรปรวนร่วมไม่สำคัญ เพราะค่าประมาณไม่เปลี่ยน"
จริงเฉพาะเมื่อข้อมูลครบถ้วนและใช้เวลาแบบจัดกลุ่ม และถึงอย่างนั้นค่าคลาดเคลื่อนมาตรฐานก็ยังต่างกัน
วิธีแก้: กำหนดแบบจำลองความแปรปรวนร่วมไว้ล่วงหน้าและรายงานคู่กับผลลัพธ์
สิ่งที่ควรทำในการวิเคราะห์ของคุณเอง
- เขียน estimand พร้อมครั้งที่วัดไว้ในโปรโตคอล
- รายงานคอนทราสต์ที่ครั้งที่วัดนั้นพร้อมช่วงเชื่อมั่น จาก
lincomหรือcontrastใน Stata หรือผลรวมของสัมประสิทธิ์แบบเดียวกันใน R - ตัดสินใจล่วงหน้าว่าจะใช้คะแนนที่จุดเริ่มต้นอย่างไร เมื่อข้อมูลครบถ้วน ANCOVA มักเป็นทางเลือกที่มีประสิทธิภาพ
- เลือกเวลาแบบจัดกลุ่ม เว้นแต่คาดว่าแนวโน้มเป็นเส้นตรง และตรวจแบบจำลองความชันกับค่าเฉลี่ยที่สังเกตได้ในแต่ละครั้งที่วัด
- ระบุแบบจำลองความแปรปรวนร่วมและข้อสมมติเรื่องข้อมูลขาดหายที่อยู่เบื้องหลัง
- ก่อนอ่านสัมประสิทธิ์ใด ตรวจว่าซอฟต์แวร์ใช้ครั้งที่วัดและกลุ่มใดเป็นระดับอ้างอิง
อภิธานศัพท์
- treatment-by-time interaction (ปฏิกิริยาสัมพันธ์ระหว่างการรักษากับเวลา)
- พจน์ในแบบจำลองที่ทำให้ความต่างระหว่างกลุ่มเปลี่ยนไปในแต่ละครั้งที่วัด
- difference in differences
- การเปลี่ยนแปลงในกลุ่มรักษาลบการเปลี่ยนแปลงในกลุ่มควบคุม
- categorical time
- แต่ละครั้งที่วัดถูกลงรหัสเป็นหมวดของตัวเอง โดยไม่สมมติรูปร่างของแนวโน้ม
- MMRM (แบบจำลองผสมสำหรับการวัดซ้ำ)
- แบบจำลองผสมสำหรับการวัดซ้ำ (mixed model for repeated measures): ครั้งที่วัดแบบจัดกลุ่ม พจน์กลุ่มการรักษาคูณครั้งที่วัด และความแปรปรวนร่วมแบบไม่มีโครงสร้างสำหรับผลลัพธ์ที่วัดซ้ำ
- unstructured covariance
- แบบจำลองความแปรปรวนร่วมที่ความแปรปรวนและสหสัมพันธ์ทุกค่าถูกประมาณอย่างอิสระ
- ANCOVA
- การวิเคราะห์ความแปรปรวนร่วม (analysis of covariance): คะแนนตอนติดตามผลที่ถดถอยบนกลุ่มและคะแนนที่จุดเริ่มต้น
- change score
- คะแนนตอนติดตามผลลบคะแนนที่จุดเริ่มต้น แล้วนำมาเปรียบเทียบระหว่างกลุ่ม
- constrained baseline
- แบบจำลองระยะยาวที่ให้ทั้งสองกลุ่มมีค่าเฉลี่ยที่จุดเริ่มต้นร่วมกันหนึ่งค่า
- estimand (ปริมาณเป้าหมายของการประมาณ)
- ปริมาณที่แน่ชัดซึ่งการศึกษาตั้งใจประมาณ ซึ่ง ICH E9(R1) อธิบายด้วยการรักษา ประชากร ตัวแปร การจัดการกับเหตุการณ์ระหว่างการศึกษา และมาตรวัดสรุประดับประชากร
เอกสารอ้างอิง
- Fitzmaurice GM, Laird NM, Ware JH. Applied longitudinal analysis. 2nd ed. Wiley; 2011. doi:10.1002/9781119513469 https://doi.org/10.1002/9781119513469
- Vickers AJ, Altman DG. Analysing controlled trials with baseline and follow up measurements. BMJ. 2001;323(7321):1123-1124. doi:10.1136/bmj.323.7321.1123 https://doi.org/10.1136/bmj.323.7321.1123
- Twisk J, Bosman L, Hoekstra T, Rijnhart J, Welten M, Heymans M. Different ways to estimate treatment effects in randomised controlled trials. Contemp Clin Trials Commun. 2018;10:80-85. doi:10.1016/j.conctc.2018.03.008 https://doi.org/10.1016/j.conctc.2018.03.008
- International Council for Harmonisation (ICH). ICH E9(R1): Addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials. ICH Harmonised Guideline, Step 4; 2019. https://database.ich.org/sites/default/files/E9-R1_Step4_Guideline_2019_1203.pdf
ประเด็นสำคัญ
- ในแบบจำลองที่มีกลุ่ม ครั้งที่วัด และปฏิกิริยาสัมพันธ์ของทั้งสอง สัมประสิทธิ์ของกลุ่มคือความต่างระหว่างกลุ่มที่ครั้งที่วัดอ้างอิง
- ความต่างระหว่างกลุ่มที่ครั้งหลังเท่ากับสัมประสิทธิ์ของกลุ่มบวกสัมประสิทธิ์ปฏิกิริยาสัมพันธ์ของครั้งนั้น
- เวลาแบบจัดกลุ่มไม่สมมติรูปร่างของแนวโน้ม ส่วนเวลาเชิงเส้นใช้ผลต่างของความชันหนึ่งค่าและสมมติว่าเป็นเส้นตรง
- ANCOVA คะแนนการเปลี่ยนแปลง และจุดเริ่มต้นร่วมแบบบังคับมุ่งไปที่ความต่างที่สัปดาห์ที่ 12 ค่าเดียวกันในการทดลองแบบสุ่ม และต่างกันเป็นหลักที่ความแม่นยำ
- ระบุครั้งที่วัดไว้ใน estimand และระบุแบบจำลองความแปรปรวนร่วมคู่กับผลลัพธ์
อ่านต่อในวิกิ: [[repeated-measures-modeling-guide-th]] [[mixed-model-series-guide-th]]