วิเคราะห์การศึกษาแบบไขว้ 2x2: แยกผลของการรักษาออกจากผลของช่วงเวลาและผลตกค้าง
On this page
Read the English version
บทคัดย่อ
การศึกษาแบบไขว้ (crossover trial) 2x2 ให้ผู้เข้าร่วมทุกคนได้รับการรักษาทั้งสองแบบ แบบละหนึ่งช่วงเวลา (period) ตามลำดับ AB หรือ BA ผลการวัดจึงขยับตามผลของการรักษา ตามผลของช่วงเวลาที่ทุกคนได้รับร่วมกันในช่วงเวลาที่ 2 และอาจขยับตามผลตกค้าง (carry-over) ซึ่งเป็นผลของยาตัวแรกที่ยังออกฤทธิ์อยู่หลังเปลี่ยนยา ครึ่งหนึ่งของความต่างระหว่างค่าเฉลี่ยของผลต่างระหว่างช่วงเวลา (period difference คือค่าในช่วงเวลาที่ 1 ลบค่าในช่วงเวลาที่ 2 ของผู้เข้าร่วมแต่ละคน) ของสองลำดับ ประมาณผลของการรักษาโดยปราศจากผลของช่วงเวลา ผลตกค้างแยกออกจากการเปรียบเทียบระหว่างลำดับไม่ได้ ในการศึกษาจำลองที่มีผู้ใหญ่ที่เป็นความดันโลหิตสูง 24 คน ค่าประมาณผลของยา A เทียบกับยา B ต่อความดันโลหิตซิสโตลิกคือ -8.02 mmHg (ช่วงเชื่อมั่น 95% คือ -11.74 ถึง -4.30) เมื่อมีผลตกค้างจำลอง -2 mmHg ตัวประมาณนี้มุ่งไปที่ -7 mmHg และการสุ่มครั้งนี้บังเอิญได้ค่าใกล้ค่าจริงที่ -8 การทดสอบผลตกค้างมีอำนาจการทดสอบ (power คือโอกาสที่จะตรวจพบผลที่มีอยู่จริง) 0.054 สิ่งที่ป้องกันผลตกค้างคือช่วงล้างยา (washout) ซึ่งเป็นช่วงเว้นการรักษาที่ยาวพอให้ฤทธิ์ของยาตัวแรกหมดไป ไม่ใช่การทดสอบเบื้องต้น
สองลำดับ สองช่วงเวลา และการเปรียบเทียบยาหนึ่งคู่
ผู้ใหญ่ที่เป็นความดันโลหิตสูงยี่สิบสี่คนเข้าร่วมการทดลองยาลดความดันโลหิตสมมติสองตัว คือยา A และยา B ทุกคนได้รับยาทั้งสองตัวสลับกัน โดยวัดความดันโลหิตซิสโตลิกเมื่อสิ้นสุดแต่ละช่วงการรักษา สิบสองคนถูกสุ่มให้ได้รับลำดับ AB และอีกสิบสองคนได้รับลำดับ BA
ในที่ประชุมคณะกรรมการกำกับการทดลอง สรุประหว่างทางแสดงว่าสองลำดับเปลี่ยนไปในทิศทางตรงข้ามกันหลังเปลี่ยนยา (ข้อมูลจำลอง) ความดันโลหิตซิสโตลิกเฉลี่ยเพิ่มขึ้น 3.88 mmHg จากช่วงเวลาที่ 1 ไปช่วงเวลาที่ 2 ในกลุ่มที่ได้ยา A ก่อน (ลำดับ AB) และลดลง 12.15 mmHg ในกลุ่มที่ได้ยา B ก่อน (ลำดับ BA) เมื่อเฉลี่ยทั้งสองลำดับ ช่วงเวลาที่ 2 อยู่ต่ำกว่า
การเปลี่ยนแปลงส่วนหนึ่งเป็นของยา ส่วนหนึ่งอาจเป็นของช่วงเวลาที่ 2 เอง และอีกส่วนอาจเป็นของสิ่งที่ยาตัวแรกทิ้งไว้ คณะกรรมการต้องการตัวเลขเดียวสำหรับผลของ A เทียบกับ B การวิเคราะห์จึงต้องแยกผลของยาออกจากช่วงเวลา และจากผลของยาตัวแรกที่ยังหลงเหลือหลังเปลี่ยนยา
การศึกษาแบบไขว้คืออะไร และอะไรบ้างที่ขยับผลการวัด
การศึกษาแบบไขว้ (crossover trial) ให้ผู้เข้าร่วมแต่ละคนได้รับทุกการรักษาที่ศึกษา ทีละแบบ ตามลำดับที่สุ่มไว้ ช่วงการรักษาแต่ละช่วงคือ ช่วงเวลา (period) และลำดับที่ผู้เข้าร่วมได้รับคือ ลำดับ (sequence) ของคนนั้น แบบแผน 2x2 มีสองช่วงเวลาและสองลำดับ คือ AB และ BA
แบบแผนนี้อาศัยการเปรียบเทียบภายในผู้เข้าร่วม (within-participant comparison) คือวัดแต่ละคนทั้งสองยา ความต่างที่คงที่ระหว่างคน เช่นระดับความดันโลหิตตามปกติ จึงหักล้างกันไป เมื่อความต่างเหล่านั้นมาก การศึกษาแบบไขว้อาจให้ความแม่นยำเท่ากับการทดลองแบบกลุ่มขนาน (parallel-group trial) ที่ใหญ่กว่ามาก ซึ่งผู้เข้าร่วมแต่ละคนได้รับการรักษาเพียงแบบเดียว [1, 2] ราคาที่ต้องจ่ายคือเวลาและลำดับของการรักษาเข้ามาอยู่ในข้อมูลด้วย
เรียกชื่อผลได้สี่อย่าง ผลของการรักษา (treatment effect) $\Delta$ (เดลตา) คือความต่างระหว่างยา A กับยา B ภายในผู้เข้าร่วมคนเดียวกัน ผลของช่วงเวลา (period effect) $\pi$ (พาย) คือการเลื่อนของค่าที่ทุกคนในช่วงเวลาที่ 2 ได้รับร่วมกัน ไม่ว่าจะได้ยาตัวใด
ผลตกค้าง (carry-over effect) $\lambda$ (แลมบ์ดา) คือผลของยาในช่วงเวลาที่ 1 ที่ยังคงอยู่ในช่วงเวลาที่ 2 ช่วงล้างยา (washout period) คือช่วงที่ไม่มีการรักษาคั่นระหว่างช่วงเวลา ยาวพอให้ผลของยาตัวแรกจางหายไป และเป็นการป้องกันผลตกค้างในขั้นออกแบบ
ผลของลำดับ (sequence effect) คือความต่างเชิงระบบใด ๆ ระหว่างกลุ่ม AB กับกลุ่ม BA ในภาพรวม ตารางแสดงว่าแต่ละผลตามอะไร และหัวข้อถัดไปแสดงเหตุผลที่สองผลหลังแยกจากกันไม่ได้
| ผล | ตามอะไร | สัญลักษณ์ | ประมาณได้หรือไม่ |
|---|---|---|---|
| ผลของการรักษา | ยาที่ได้รับในช่วงเวลาปัจจุบัน | $\Delta$ | ได้ ภายในผู้เข้าร่วม (เลื่อนไปลบครึ่งหนึ่งของผลตกค้าง ถ้ามี) |
| ผลของช่วงเวลา | ค่าที่วัดอยู่ในช่วงเวลาที่ 1 หรือช่วงเวลาที่ 2 | $\pi$ | ได้ ภายในผู้เข้าร่วม (เลื่อนไปบวกครึ่งหนึ่งของผลตกค้าง ถ้ามี) |
| ผลตกค้าง | ยาที่ได้รับในช่วงเวลาก่อนหน้า | $\lambda$ | ได้เฉพาะระหว่างผู้เข้าร่วม และได้อย่างไม่แม่นยำ |
| ผลของลำดับ | ลำดับ AB หรือ BA | ไม่มีสัญลักษณ์ของตัวเอง | ไม่ได้ แบบแผนทำให้มันมีรูปแบบเดียวกันในข้อมูลกับผลตกค้าง (ทั้งสองปนกัน หรือ aliased) จึงไม่มีการวิเคราะห์ใดแยกได้ |
แบบจำลองสำหรับค่าเฉลี่ยของเซลล์ทั้งสี่
ให้ $Y_{ij}$ เป็นผลลัพธ์ของผู้เข้าร่วมคนที่ $i$ ในช่วงเวลาที่ $j$ โดย $j$ เป็น 1 หรือ 2 แบบแผนนี้มีสี่เซลล์ คือหนึ่งเซลล์ต่อหนึ่งลำดับในหนึ่งช่วงเวลา และค่าเฉลี่ยของเซลล์ (cell mean) คือค่าเฉลี่ยของการวัดในเซลล์หนึ่ง แบบจำลองมาตรฐานของการศึกษาแบบไขว้ 2x2 คือ
$$Y_{ij} = \mu + s_i + \pi\,[\text{period 2}] + \Delta\,[\text{on A}] + \lambda\,[\text{period 2 of AB}] + e_{ij}$$วงเล็บเหลี่ยมแต่ละอันมีค่า 1 เมื่อเงื่อนไขในวงเล็บเป็นจริง และ 0 ในกรณีอื่น ค่าคงที่ $\mu$ (มิว) คือค่าเฉลี่ยของยา B ในช่วงเวลาที่ 1 $s_i$ คือระดับของผู้เข้าร่วมคนที่ $i$ เมื่อเทียบกับค่านั้น และ $e_{ij}$ คือความคลาดเคลื่อนของการวัดครั้งเดียว การมอง $s_i$ เป็นค่าที่สุ่มมาจากการแจกแจงปกติทำให้มันเป็นจุดตัดแกนสุ่ม (random intercept) ซึ่งเป็นพจน์ที่แบบจำลองผสม (mixed model คือการถดถอยที่มีระดับสุ่มของผู้เข้าร่วมแต่ละคน) ใช้แทนระดับของผู้เข้าร่วมแต่ละคน
พูดให้เคร่งครัด $\lambda$ คือผลตกค้างของ A ลบผลตกค้างของ B ผลตกค้างที่เกิดเท่ากันกับยาทั้งสองตัวดูเหมือนผลของช่วงเวลาทุกประการ แบบจำลองจึงมีเพียงส่วนต่าง โดยวางไว้ที่ช่วงเวลาที่ 2 ของลำดับ AB
| ลำดับ | ช่วงเวลาที่ 1 | ช่วงเวลาที่ 2 |
|---|---|---|
| AB | $\mu + \Delta$ (ได้ยา A) | $\mu + \pi + \lambda$ (ได้ยา B) |
| BA | $\mu$ (ได้ยา B) | $\mu + \pi + \Delta$ (ได้ยา A) |
ผลต่างระหว่างช่วงเวลา: ตัวประมาณ
สำหรับผู้เข้าร่วมแต่ละคน ให้นำค่าในช่วงเวลาที่ 1 ลบด้วยค่าในช่วงเวลาที่ 2 คือ $d_i = Y_{i1} - Y_{i2}$ ระดับ $s_i$ ของผู้เข้าร่วมปรากฏในการวัดทั้งสองครั้งและหักล้างกันไป ให้ $\bar d_{AB}$ และ $\bar d_{BA}$ เป็นค่าเฉลี่ยของ $d_i$ ในแต่ละลำดับ
จากค่าเฉลี่ยของเซลล์ $\bar d_{AB}$ ประมาณ $\Delta - \pi - \lambda$ และ $\bar d_{BA}$ ประมาณ $-\Delta - \pi$ ผลต่างของสองค่านี้ตัด $\pi$ ออก และผลรวมของมันตัด $\Delta$ ออก:
$$\hat\Delta = \tfrac{1}{2}\left(\bar d_{AB} - \bar d_{BA}\right)$$ตัวประมาณผลของการรักษา $\hat\Delta$ (เดลตาแฮต เครื่องหมายหมวกแสดงว่าเป็นค่าประมาณ) คือครึ่งหนึ่งของผลต่างระหว่างค่าเฉลี่ยของผลต่างระหว่างช่วงเวลาของสองลำดับ
$$\hat\pi = -\tfrac{1}{2}\left(\bar d_{AB} + \bar d_{BA}\right)$$ตัวประมาณผลของช่วงเวลา $\hat\pi$ คือลบครึ่งหนึ่งของผลรวมนั้น ซึ่งเป็นการเลื่อนจากช่วงเวลาที่ 1 ไปช่วงเวลาที่ 2 เมื่อไม่มีผลตกค้าง ($\lambda = 0$) $\hat\Delta$ ประมาณ $\Delta$ และ $\hat\pi$ ประมาณ $\pi$
ค่าคลาดเคลื่อนมาตรฐาน (standard error, SE) ของ $\hat\Delta$ เท่ากับครึ่งหนึ่งของค่าคลาดเคลื่อนมาตรฐานของการทดสอบที (t-test) ของตัวอย่างสองกลุ่มที่เปรียบเทียบ $d_i$ ระหว่างลำดับ (Hills และ Armitage [3]) ค่าคลาดเคลื่อนมาตรฐานคือการกระจายของค่าประมาณเมื่อทำการทดลองเดิมซ้ำหลายครั้ง การทดสอบทีนั้นประมาณ $\bar d_{AB} - \bar d_{BA}$ ซึ่งเป็นสองเท่าของ $\hat\Delta$ ค่าสถิติทีและค่า P ของมันจึงใช้กับ $\hat\Delta$ ได้โดยไม่ต้องเปลี่ยน
ตัวอย่างคำนวณด้วยมือ: มีผลของช่วงเวลา แต่ยาไม่ต่างกัน
ตัวอย่างคำนวณด้วยมือนี้ใช้คะแนนอาการตั้งแต่ 0 ถึง 10 โดยคะแนนสูงหมายถึงแย่กว่า ในลำดับ AB คะแนนเฉลี่ยเท่ากับ 8 เมื่อได้ยา A ในช่วงเวลาที่ 1 และ 5 เมื่อได้ยา B ในช่วงเวลาที่ 2 ในลำดับ BA คะแนนเฉลี่ยเท่ากับ 8 เมื่อได้ยา B ในช่วงเวลาที่ 1 และ 5 เมื่อได้ยา A ในช่วงเวลาที่ 2
-
ผลต่างระหว่างช่วงเวลาเฉลี่ยในลำดับ AB
\[ \bar d_{AB} = 8 - 5 = 3 \]
ช่วงเวลาที่ 1 ลบช่วงเวลาที่ 2 ในลำดับ AB
-
ผลต่างระหว่างช่วงเวลาเฉลี่ยในลำดับ BA
\[ \bar d_{BA} = 8 - 5 = 3 \]
การลบแบบเดียวกันในลำดับ BA
-
ผลของการรักษา
\[ \hat\Delta = \tfrac{1}{2}(3 - 3) = 0 \]
ยา A กับยา B ไม่ต่างกัน
-
ผลของช่วงเวลา
\[ \hat\pi = -\tfrac{1}{2}(3 + 3) = -3 \]
คะแนนในช่วงเวลาที่ 2 ต่ำกว่า 3 คะแนน ไม่ว่าจะได้ยาตัวใด
ผลลัพธ์: การลดลงทั้งหมดจาก 8 เป็น 5 เป็นของช่วงเวลาที่ 2 และไม่ได้เป็นของยาตัวใดเลย
เมื่อดูเฉพาะภายในลำดับ AB การลดลงเกิดพร้อมกับการเปลี่ยนไปใช้ยา B และจะถูกยกให้เป็นผลของ B เมื่อดูเฉพาะภายในลำดับ BA การลดลงเดียวกันจะถูกยกให้เป็นผลของ A มีเพียงการนำสองลำดับมาพิจารณาร่วมกันเท่านั้นที่แยกยาออกจากช่วงเวลาได้
ผลตกค้างปนอยู่กับลำดับ
ทีนี้ให้ยาตัวแรกยังหลงเหลืออยู่ นำค่าเฉลี่ยของเซลล์ใส่ในตัวประมาณทั้งสองแล้วหาค่าคาดหมาย เขียนเป็น $E(\cdot)$ คือค่าเฉลี่ยเมื่อทำการทดลองเดิมซ้ำหลายครั้ง:
$$E(\hat\Delta) = \Delta - \tfrac{\lambda}{2}$$ตัวประมาณผลของการรักษาจึงเลื่อนไปลบครึ่งหนึ่งของผลตกค้าง ตัวประมาณผลของช่วงเวลาเลื่อนไปอีกทาง คือ $E(\hat\pi) = \pi + \lambda/2$
ผลต่างระหว่างช่วงเวลาเฉลี่ยทั้งสองค่าเป็นข้อมูลภายในผู้เข้าร่วมเพียงอย่างเดียว และให้สองสมการ คือ $\Delta - \pi - \lambda$ และ $-\Delta - \pi$ สำหรับสามตัวไม่ทราบค่า ไม่มีการรวมใดของมันที่แยก $\Delta$ หรือ $\pi$ ออกจาก $\lambda$ ได้ ตัวประมาณเดียวของ $\Delta$ ที่ปราศจาก $\lambda$ คือการเปรียบเทียบช่วงเวลาที่ 1 อย่างเดียว (AB ช่วงเวลาที่ 1 ลบ BA ช่วงเวลาที่ 1) ซึ่งเป็นการเปรียบเทียบระหว่างผู้เข้าร่วม แบบจำลองที่มีพจน์ผลตกค้างในหัวข้อหลังของบทความนี้ทำเช่นนั้นพอดี
สองผลถือว่าปนกัน (aliased) เมื่อแบบแผนทำให้ทั้งสองมีรูปแบบเดียวกันในข้อมูล จึงไม่มีการวิเคราะห์ใดแยกได้
ในการศึกษาแบบไขว้ 2x2 ผลตกค้าง ผลของลำดับ และปฏิกิริยาสัมพันธ์ระหว่างการรักษากับช่วงเวลา (treatment-by-period interaction คือความต่างของยาที่เปลี่ยนไประหว่างช่วงเวลา) ล้วนปรากฏเป็นปริมาณเดียวกัน คือความต่างระหว่างลำดับของผลรวมรายผู้เข้าร่วม (ค่าช่วงเวลาที่ 1 บวกค่าช่วงเวลาที่ 2 ของแต่ละคน) [1, 2] ข้อมูลบอกไม่ได้ว่าปริมาณนี้วัดผลตัวใดในสามตัวนี้
การศึกษาจำลองด้านล่างสร้างด้วย $\Delta = -8$, $\pi = -2$ และ $\lambda = -2$ mmHg ตัวประมาณผลของการรักษาจึงมุ่งไปที่ $-8 - (-2)/2 = -7$ mmHg ซึ่งเอนเอียง 1 mmHg เข้าหาการไม่มีผล ส่วนตัวประมาณผลของช่วงเวลามุ่งไปที่ $-2 + (-2)/2 = -3$ mmHg
ทำไมการทดสอบผลตกค้างจึงมีอำนาจการทดสอบต่ำ
ผลตกค้างประมาณได้จากผลรวมรายผู้เข้าร่วมเท่านั้น คือ $t_i = Y_{i1} + Y_{i2}$ ค่าเฉลี่ยของมันต่างกันระหว่างลำดับด้วย $\lambda$ และการทดสอบทีของตัวอย่างสองกลุ่มบนผลรวมคือการทดสอบผลตกค้าง อย่างไรก็ตาม ผลรวมแต่ละค่ามี $2s_i$ อยู่ด้วย ความแปรผันระหว่างผู้เข้าร่วมทั้งหมดจึงยังคงอยู่ในนั้น:
$$\operatorname{SE}(\hat\lambda) = \sqrt{\left(4\sigma_s^2 + 2\sigma_e^2\right)\left(\tfrac{1}{n_{AB}} + \tfrac{1}{n_{BA}}\right)}$$ในสมการนี้ $\sigma_s$ คือส่วนเบี่ยงเบนมาตรฐาน (SD) ระหว่างผู้เข้าร่วม ซึ่งเป็นการกระจายของระดับ $s_i$ และ $\sigma_e$ คือ SD ภายในผู้เข้าร่วมของการวัดครั้งเดียว ส่วน $n_{AB}$ และ $n_{BA}$ คือจำนวนผู้เข้าร่วมในแต่ละลำดับ ตัวประมาณผลของการรักษาทำงานบนผลต่าง ซึ่ง $\sigma_s$ หักล้างกัน ค่าคลาดเคลื่อนมาตรฐานของมันจึงขึ้นกับ $\sigma_e$ เท่านั้น
ด้วยค่าที่ใช้ออกแบบการศึกษาจำลอง คือ $\sigma_s = 12$ mmHg, $\sigma_e = 6$ mmHg และผู้เข้าร่วม 12 คนต่อลำดับ ค่าคลาดเคลื่อนมาตรฐานของตัวประมาณผลตกค้างคือ 10.39 mmHg
อำนาจการทดสอบ (power) คือความน่าจะเป็นที่การทดสอบจะประกาศว่ามีผล เมื่อมีผลขนาดที่ระบุอยู่จริง อำนาจการทดสอบของการทดสอบสองทางที่ระดับ 5% ต่อผลตกค้างจริง $\lambda = -2$ mmHg คือ 0.054 ซึ่งสูงกว่า 5% ที่จะได้เมื่อไม่มีผลตกค้างเลยเพียงนิดเดียว ผลตกค้างที่แบบแผนนี้จะตรวจพบได้ในสี่จากห้าการศึกษา ซึ่งเป็นเป้าหมายอำนาจการทดสอบที่มักใช้ในการวางแผน คือประมาณ 30.5 mmHg ซึ่งใหญ่กว่าผลของการรักษามาก
การทดสอบผลตกค้างเปรียบเทียบผลรวมของผู้ร่วมวิจัยแต่ละคนระหว่างลำดับการได้รับยา ซึ่งเป็นการเปรียบเทียบระหว่างบุคคลที่มี power ต่ำ ผลที่ไม่มีนัยสำคัญจึงไม่ได้แสดงว่าไม่มีผลตกค้าง ผลตกค้างป้องกันได้ด้วยช่วงล้างยาที่เพียงพอ ไม่ใช่กำจัดได้ด้วยการทดสอบเบื้องต้น
ทดสอบก่อนแล้วใช้เฉพาะช่วงเวลาที่ 1: ทำไมขั้นตอนสองขั้นจึงล้มเหลว
แนวทางที่เคยเป็นมาตรฐาน ซึ่ง Grizzle เสนอ [4] มีสองขั้น ขั้นแรกทดสอบผลตกค้างบนผลรวมรายผู้เข้าร่วมที่ระดับนัยสำคัญที่ผ่อนปรน ถ้าการทดสอบนั้นมีนัยสำคัญ ให้ทิ้งช่วงเวลาที่ 2 แล้วเปรียบเทียบยาในช่วงเวลาที่ 1 เหมือนการทดลองแบบกลุ่มขนาน ถ้าไม่เช่นนั้นให้ใช้ผลต่างระหว่างช่วงเวลา
Freeman [5] แสดงว่าขั้นตอนนี้ทำให้ความคลาดเคลื่อนประเภทที่ 1 (type I error) ซึ่งคือความน่าจะเป็นที่จะประกาศว่ายาต่างกันทั้งที่ไม่ต่างกันเลย สูงเกินจริง และทำให้ค่าประมาณสุดท้ายเอนเอียง การทดสอบผลตกค้างกับการเปรียบเทียบช่วงเวลาที่ 1 ใช้ข้อมูลช่วงเวลาที่ 1 ร่วมกันและสัมพันธ์กันมาก ขั้นแรกจึงเป็นตัวกำหนดว่าขั้นที่สองจะรายงานผลแบบไหน
งานด้านระเบียบวิธีในปัจจุบันแนะนำให้เลี่ยงขั้นตอนสองขั้น [1, 5] ผลตกค้างถูกจัดการในขั้นออกแบบแทน ด้วยภาวะที่คงที่ตลอดการทดลอง ช่วงล้างยาระหว่างการรักษาที่ยาวพอให้ผลของยาตัวแรกจางหาย และการวัดผลลัพธ์เมื่อสิ้นสุดแต่ละช่วงเวลา คู่มือการออกแบบและวิเคราะห์การศึกษาแบบไขว้ (หน้าภาษาอังกฤษ) ครอบคลุมด้านการออกแบบ รวมถึงช่วงล้างยา อย่างลึกกว่านี้
แบบจำลองสำหรับวิเคราะห์: แบบจำลองผสมที่มีการรักษาและช่วงเวลา
การวิเคราะห์ที่ใช้กันทั่วไปคือแบบจำลองผสม (mixed model) เชิงเส้น ซึ่งเป็นการถดถอยที่มีพจน์คงที่สำหรับการรักษาและช่วงเวลา และจุดตัดแกนสุ่มสำหรับผู้เข้าร่วมแต่ละคน ประมาณด้วย REML (restricted maximum likelihood หรือภาวะน่าจะเป็นสูงสุดแบบจำกัด ซึ่งเป็นวิธีที่ใช้กันทั่วไปในการประมาณการกระจายระหว่างผู้เข้าร่วมและภายในผู้เข้าร่วมในตัวอย่างขนาดเล็ก)
องศาอิสระ (degrees of freedom) ของการทดสอบที กำหนดรูปร่างที่แน่นอนของการแจกแจงทีที่อยู่เบื้องหลังค่า P และช่วงเชื่อมั่น การทดสอบทีของแบบจำลองนี้ใช้องศาอิสระแบบ Kenward-Roger ซึ่งเป็นการปรับสำหรับตัวอย่างขนาดเล็ก
เมื่อข้อมูลครบถ้วน สัมประสิทธิ์การรักษาของแบบจำลองนี้เท่ากับ $\hat\Delta$ พอดี และค่าคลาดเคลื่อนมาตรฐานของมันเท่ากับครึ่งหนึ่งของค่าคลาดเคลื่อนมาตรฐานของการทดสอบทีสองกลุ่มแบบความแปรปรวนรวม (pooled-variance) หรือการทดสอบทีแบบ Student บนผลต่างระหว่างช่วงเวลา ซึ่งเป็นแบบที่สมมติว่าสองลำดับมีความแปรปรวนร่วมกันหนึ่งค่า ค่าสถิติที ค่า P และองศาอิสระ 22 เหมือนกันทุกประการ การทดสอบทีแบบ Welch ซึ่งไม่รวมความแปรปรวนของสองกลุ่มจะให้องศาอิสระต่างไปเล็กน้อย
แบบจำลองแสดงคุณค่าเมื่อข้อมูลไม่เรียบร้อย มันรับตัวแปรร่วมที่จุดเริ่มต้น (baseline covariates) ได้ และเก็บผู้เข้าร่วมที่ขาดข้อมูลหนึ่งช่วงเวลาไว้ได้ ชุดบทความแบบจำลองผสม ครอบคลุมแบบจำลองนี้อย่างลึก
การเพิ่มพจน์ผลตกค้างไม่ได้ช่วยกู้การวิเคราะห์ เมื่อมี $\lambda$ ในแบบจำลอง สัมประสิทธิ์การรักษาจะถูกประมาณจากช่วงเวลาที่ 1 อย่างเดียว ซึ่งเป็นการเปรียบเทียบระหว่างผู้เข้าร่วมที่สละข้อได้เปรียบหลักของการศึกษาแบบไขว้
การศึกษาจำลองใน Stata และ R
ชุดข้อมูล (ข้อมูลจำลอง) มี 48 แถว หนึ่งแถวต่อผู้เข้าร่วมหนึ่งคนต่อหนึ่งช่วงเวลา คือผู้เข้าร่วม 24 คน ลำดับละ 12 คน โดยวัดความดันโลหิตซิสโตลิกเป็น mmHg สร้างด้วย $\Delta = -8$, $\pi = -2$ และ $\lambda = -2$ mmHg ส่วนเบี่ยงเบนมาตรฐานระหว่างผู้เข้าร่วม 12 mmHg และส่วนเบี่ยงเบนมาตรฐานภายในผู้เข้าร่วม 6 mmHg
ผลต่างระหว่างช่วงเวลาเฉลี่ยคือ -3.88 mmHg ในลำดับ AB และ 12.15 mmHg ในลำดับ BA ครึ่งหนึ่งของผลต่างของสองค่านี้ให้ $\hat\Delta = -8.02$ mmHg (ค่าคลาดเคลื่อนมาตรฐาน 1.79 ช่วงเชื่อมั่น 95% คือ -11.74 ถึง -4.30 ค่า P = 0.0002) แบบจำลองผสมให้สัมประสิทธิ์ ค่าคลาดเคลื่อนมาตรฐาน และช่วงเชื่อมั่นเดียวกัน ค่าประมาณผลของช่วงเวลาคือ -4.13 mmHg (ช่วงเชื่อมั่น 95% คือ -7.85 ถึง -0.41)
ในการสุ่มครั้งเดียวนี้ $\hat\Delta$ บังเอิญอยู่ใกล้ $\Delta$ จริงที่ -8 แม้เป้าหมายของมันคือ -7 การศึกษาเดียวที่มีผู้เข้าร่วม 24 คนไม่อาจเผยความเอนเอียง 1 mmHg เมื่อเทียบกับค่าคลาดเคลื่อนมาตรฐาน 1.79 ได้
ค่าประมาณผลตกค้างคือ 7.45 mmHg (ค่าคลาดเคลื่อนมาตรฐาน 5.45 ช่วงเชื่อมั่น 95% คือ -3.85 ถึง 18.75 ค่า P = 0.19) เทียบกับค่าจริงที่ -2 mmHg ซึ่งห่างจากความจริงมากและไม่มีนัยสำคัญ ตรงกับที่อำนาจการทดสอบต่ำทำนายไว้ ค่าคลาดเคลื่อนมาตรฐานของมันมีขนาดราวครึ่งหนึ่งของค่าตามการออกแบบที่ 10.39 เพราะส่วนเบี่ยงเบนมาตรฐานระหว่างผู้เข้าร่วมที่ประมาณได้ในการสุ่มครั้งนี้ คือ 5.19 mmHg ต่ำกว่า 12 mmHg ที่ใช้สร้างข้อมูลมาก
แบบจำลองที่มีพจน์ผลตกค้างประมาณผลของการรักษาจากช่วงเวลาที่ 1 อย่างเดียว ได้ -4.29 mmHg (ช่วงเชื่อมั่น 95% คือ -10.89 ถึง 2.31)
| ปริมาณ | ค่าประมาณ (mmHg) | ช่วงเชื่อมั่น 95% (mmHg) | เป้าหมาย (mmHg) |
|---|---|---|---|
| ผลต่างระหว่างช่วงเวลาเฉลี่ย ลำดับ AB | -3.88 | ไม่แสดง | -4 |
| ผลต่างระหว่างช่วงเวลาเฉลี่ย ลำดับ BA | 12.15 | ไม่แสดง | 10 |
| ผลของการรักษา จากผลต่างระหว่างช่วงเวลา | -8.02 | -11.74 ถึง -4.30 | -7 |
| ผลของการรักษา จากแบบจำลองผสม | -8.02 | -11.74 ถึง -4.30 | -7 |
| ผลของช่วงเวลา | -4.13 | -7.85 ถึง -0.41 | -3 |
| ผลตกค้าง จากผลรวม AB ลบ BA | 7.45 | -3.85 ถึง 18.75 | -2 |
| ผลของการรักษา จากแบบจำลองที่มีพจน์ผลตกค้าง | -4.29 | -10.89 ถึง 2.31 | -8 |
Stata: การวิเคราะห์ทั้งหมด ตั้งแต่ค่าเฉลี่ยของเซลล์จนถึงอำนาจการทดสอบของการทดสอบผลตกค้าง
* Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any real drug)
* Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 12 participants, two periods.
version 18
clear all
set more off
set linesize 120
* read the simulated dataset (not published with this article)
import delimited using crossover.csv, clear varnames(1)
describe, short
list in 1/6, noobs sepby(id)
* participants per sequence (one row per participant)
egen byte first = tag(id)
tabulate seq if first
* cell means by sequence and period
table seq period, statistic(mean sbp) nformat(%9.2f)
* mixed model: treatment and period as fixed effects, a random intercept per participant
* (REML, Kenward-Roger degrees of freedom). Fitted quietly; the lines below print its treatment
* and period rows and the two standard deviations. R prints the same model in full.
quietly mixed sbp i.trt i.period || id:, reml dfmethod(kroger)
matrix T = r(table)
scalar mm_d = T[rownumb(T, "b"), colnumb(T, "sbp:1.trt")]
scalar mm_se = T[rownumb(T, "se"), colnumb(T, "sbp:1.trt")]
scalar mm_lo = T[rownumb(T, "ll"), colnumb(T, "sbp:1.trt")]
scalar mm_hi = T[rownumb(T, "ul"), colnumb(T, "sbp:1.trt")]
scalar mm_df = T[rownumb(T, "df"), colnumb(T, "sbp:1.trt")]
scalar mm_p = T[rownumb(T, "pvalue"), colnumb(T, "sbp:1.trt")]
scalar mm_pi = T[rownumb(T, "b"), colnumb(T, "sbp:2.period")]
scalar mm_pilo = T[rownumb(T, "ll"), colnumb(T, "sbp:2.period")]
scalar mm_pihi = T[rownumb(T, "ul"), colnumb(T, "sbp:2.period")]
scalar mm_sdid = exp(_b[lns1_1_1:_cons])
scalar mm_sdres = exp(_b[lnsig_e:_cons])
display "mixed model, treatment (A minus B): " %8.4f mm_d " SE " %7.4f mm_se " 95% CI " %8.4f mm_lo " to " %8.4f mm_hi
display "mixed model, treatment: df " %7.4f mm_df " P " %6.4f mm_p
display "mixed model, period (2 minus 1): " %8.4f mm_pi " 95% CI " %8.4f mm_pilo " to " %8.4f mm_pihi
display "mixed model, SD between participants " %7.4f mm_sdid " SD within participants " %7.4f mm_sdres
* mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
quietly mixed sbp i.trt i.period i.carry || id:, reml dfmethod(kroger)
matrix C = r(table)
scalar mmc_d = C[rownumb(C, "b"), colnumb(C, "sbp:1.trt")]
scalar mmc_lo = C[rownumb(C, "ll"), colnumb(C, "sbp:1.trt")]
scalar mmc_hi = C[rownumb(C, "ul"), colnumb(C, "sbp:1.trt")]
scalar mmc_l = C[rownumb(C, "b"), colnumb(C, "sbp:1.carry")]
scalar mmc_lse = C[rownumb(C, "se"), colnumb(C, "sbp:1.carry")]
scalar mmc_lp = C[rownumb(C, "pvalue"), colnumb(C, "sbp:1.carry")]
display "carry-over model, treatment (A minus B): " %8.4f mmc_d " 95% CI " %8.4f mmc_lo " to " %8.4f mmc_hi
display "carry-over model, carry-over term: " %8.4f mmc_l " SE " %7.4f mmc_lse " P " %6.4f mmc_lp
* one row per participant: period difference d = period 1 minus period 2, and the participant total
keep id seq period sbp
reshape wide sbp, i(id) j(period)
generate double d = sbp1 - sbp2
generate double tot = sbp1 + sbp2
generate byte grp = cond(seq == "AB", 1, 2)
label define grp 1 "AB" 2 "BA"
label values grp grp
* treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
ttest d, by(grp)
scalar dbar_ab = r(mu_1)
scalar dbar_ba = r(mu_2)
scalar df_d = r(df_t)
scalar delta = (r(mu_1) - r(mu_2)) / 2
scalar delta_se = r(se) / 2
scalar delta_p = r(p)
scalar tc = invttail(df_d, 0.025)
display "treatment (A minus B): " %8.4f delta " SE " %7.4f delta_se " 95% CI " %8.4f delta - tc * delta_se " to " %8.4f delta + tc * delta_se " P " %6.4f delta_p
* period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
scalar per = -(dbar_ab + dbar_ba) / 2
display "period (2 minus 1): " %8.4f per " 95% CI " %8.4f per - tc * delta_se " to " %8.4f per + tc * delta_se
* carry-over test: compare participant totals between sequences (between-participant, so low power)
ttest tot, by(grp)
* design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
* power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
* (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
scalar sig_s = 12
scalar sig_e = 6
scalar lam = -2
count if grp == 1
scalar n_ab = r(N)
count if grp == 2
scalar n_ba = r(N)
scalar des_se = sqrt((4 * sig_s^2 + 2 * sig_e^2) * (1 / n_ab + 1 / n_ba))
scalar des_df = n_ab + n_ba - 2
scalar des_t = invttail(des_df, 0.025)
scalar des_power = nt(des_df, lam / des_se, -des_t) + nttail(des_df, lam / des_se, des_t)
scalar des_mde80 = (des_t + invttail(des_df, 0.20)) * des_se
display "design SE of the carry-over estimate: " %7.4f des_se
display "power of the carry-over test at lambda = -2: " %6.4f des_power
display "carry-over detected with 80% power: " %7.4f des_mde80. * Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any r
> eal drug)
. * Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 1
> 2 participants, two periods.
. version 18
. clear all
. set more off
. set linesize 120
.
. * read the simulated dataset (not published with this article)
. import delimited using crossover.csv, clear varnames(1)
(encoding automatically selected: ISO-8859-1)
(7 vars, 48 obs)
. describe, short
Contains data
Observations: 48
Variables: 7
Sorted by:
Note: Dataset has changed since last saved.
. list in 1/6, noobs sepby(id)
+----------------------------------------------------+
| id seq period treatm~t trt carry sbp |
|----------------------------------------------------|
| 1 BA 1 B 0 0 146.2 |
| 1 BA 2 A 1 0 134.5 |
|----------------------------------------------------|
| 2 AB 1 A 1 0 152.4 |
| 2 AB 2 B 0 1 149.6 |
|----------------------------------------------------|
| 3 AB 1 A 1 0 160.7 |
| 3 AB 2 B 0 1 156.8 |
+----------------------------------------------------+
.
. * participants per sequence (one row per participant)
. egen byte first = tag(id)
. tabulate seq if first
seq | Freq. Percent Cum.
------------+-----------------------------------
AB | 12 50.00 50.00
BA | 12 50.00 100.00
------------+-----------------------------------
Total | 24 100.00
.
. * cell means by sequence and period
. table seq period, statistic(mean sbp) nformat(%9.2f)
-----------------------------------
| period
| 1 2 Total
--------+--------------------------
seq |
AB | 145.72 149.60 147.66
BA | 150.01 137.86 143.93
Total | 147.86 143.73 145.80
-----------------------------------
.
. * mixed model: treatment and period as fixed effects, a random intercept per participant
. * (REML, Kenward-Roger degrees of freedom). Fitted quietly; the lines below print its treatment
. * and period rows and the two standard deviations. R prints the same model in full.
. quietly mixed sbp i.trt i.period || id:, reml dfmethod(kroger)
. matrix T = r(table)
. scalar mm_d = T[rownumb(T, "b"), colnumb(T, "sbp:1.trt")]
. scalar mm_se = T[rownumb(T, "se"), colnumb(T, "sbp:1.trt")]
. scalar mm_lo = T[rownumb(T, "ll"), colnumb(T, "sbp:1.trt")]
. scalar mm_hi = T[rownumb(T, "ul"), colnumb(T, "sbp:1.trt")]
. scalar mm_df = T[rownumb(T, "df"), colnumb(T, "sbp:1.trt")]
. scalar mm_p = T[rownumb(T, "pvalue"), colnumb(T, "sbp:1.trt")]
. scalar mm_pi = T[rownumb(T, "b"), colnumb(T, "sbp:2.period")]
. scalar mm_pilo = T[rownumb(T, "ll"), colnumb(T, "sbp:2.period")]
. scalar mm_pihi = T[rownumb(T, "ul"), colnumb(T, "sbp:2.period")]
. scalar mm_sdid = exp(_b[lns1_1_1:_cons])
. scalar mm_sdres = exp(_b[lnsig_e:_cons])
. display "mixed model, treatment (A minus B): " %8.4f mm_d " SE " %7.4f mm_se " 95% CI " %8.4f mm_lo " to " %8.4f mm_
> hi
mixed model, treatment (A minus B): -8.0167 SE 1.7938 95% CI -11.7367 to -4.2966
. display "mixed model, treatment: df " %7.4f mm_df " P " %6.4f mm_p
mixed model, treatment: df 22.0000 P 0.0002
. display "mixed model, period (2 minus 1): " %8.4f mm_pi " 95% CI " %8.4f mm_pilo " to " %8.4f mm_pihi
mixed model, period (2 minus 1): -4.1333 95% CI -7.8534 to -0.4133
. display "mixed model, SD between participants " %7.4f mm_sdid " SD within participants " %7.4f mm_sdres
mixed model, SD between participants 5.1869 SD within participants 6.2139
.
. * mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
. quietly mixed sbp i.trt i.period i.carry || id:, reml dfmethod(kroger)
. matrix C = r(table)
. scalar mmc_d = C[rownumb(C, "b"), colnumb(C, "sbp:1.trt")]
. scalar mmc_lo = C[rownumb(C, "ll"), colnumb(C, "sbp:1.trt")]
. scalar mmc_hi = C[rownumb(C, "ul"), colnumb(C, "sbp:1.trt")]
. scalar mmc_l = C[rownumb(C, "b"), colnumb(C, "sbp:1.carry")]
. scalar mmc_lse = C[rownumb(C, "se"), colnumb(C, "sbp:1.carry")]
. scalar mmc_lp = C[rownumb(C, "pvalue"), colnumb(C, "sbp:1.carry")]
. display "carry-over model, treatment (A minus B): " %8.4f mmc_d " 95% CI " %8.4f mmc_lo " to " %8.4f mmc_hi
carry-over model, treatment (A minus B): -4.2917 95% CI -10.8943 to 2.3110
. display "carry-over model, carry-over term: " %8.4f mmc_l " SE " %7.4f mmc_lse " P " %6.4f mmc_lp
carry-over model, carry-over term: 7.4500 SE 5.4483 P 0.1853
.
. * one row per participant: period difference d = period 1 minus period 2, and the participant total
. keep id seq period sbp
. reshape wide sbp, i(id) j(period)
(j = 1 2)
Data Long -> Wide
-----------------------------------------------------------------------------
Number of observations 48 -> 24
Number of variables 4 -> 4
j variable (2 values) period -> (dropped)
xij variables:
sbp -> sbp1 sbp2
-----------------------------------------------------------------------------
. generate double d = sbp1 - sbp2
. generate double tot = sbp1 + sbp2
. generate byte grp = cond(seq == "AB", 1, 2)
. label define grp 1 "AB" 2 "BA"
. label values grp grp
.
. * treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
. ttest d, by(grp)
Two-sample t test with equal variances
------------------------------------------------------------------------------
Group | Obs Mean Std. err. Std. dev. [95% conf. interval]
---------+--------------------------------------------------------------------
AB | 12 -3.883336 2.973998 10.30223 -10.42906 2.66239
BA | 12 12.15 2.006485 6.950669 7.733753 16.56624
---------+--------------------------------------------------------------------
Combined | 24 4.133331 2.423217 11.87129 -.8794746 9.146137
---------+--------------------------------------------------------------------
diff | -16.03333 3.587569 -23.4735 -8.593171
------------------------------------------------------------------------------
diff = mean(AB) - mean(BA) t = -4.4691
H0: diff = 0 Degrees of freedom = 22
Ha: diff < 0 Ha: diff != 0 Ha: diff > 0
Pr(T < t) = 0.0001 Pr(|T| > |t|) = 0.0002 Pr(T > t) = 0.9999
. scalar dbar_ab = r(mu_1)
. scalar dbar_ba = r(mu_2)
. scalar df_d = r(df_t)
. scalar delta = (r(mu_1) - r(mu_2)) / 2
. scalar delta_se = r(se) / 2
. scalar delta_p = r(p)
. scalar tc = invttail(df_d, 0.025)
. display "treatment (A minus B): " %8.4f delta " SE " %7.4f delta_se " 95% CI " %8.4f delta - tc * delta_se " to " %8
> .4f delta + tc * delta_se " P " %6.4f delta_p
treatment (A minus B): -8.0167 SE 1.7938 95% CI -11.7367 to -4.2966 P 0.0002
.
. * period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
. scalar per = -(dbar_ab + dbar_ba) / 2
. display "period (2 minus 1): " %8.4f per " 95% CI " %8.4f per - tc * delta_se " to " %8.4f per + tc * delta_se
period (2 minus 1): -4.1333 95% CI -7.8534 to -0.4132
.
. * carry-over test: compare participant totals between sequences (between-participant, so low power)
. ttest tot, by(grp)
Two-sample t test with equal variances
------------------------------------------------------------------------------
Group | Obs Mean Std. err. Std. dev. [95% conf. interval]
---------+--------------------------------------------------------------------
AB | 12 295.3167 3.709261 12.84926 287.1526 303.4807
BA | 12 287.8667 3.990715 13.82424 279.0832 296.6502
---------+--------------------------------------------------------------------
Combined | 24 291.5917 2.7752 13.59565 285.8507 297.3326
---------+--------------------------------------------------------------------
diff | 7.449999 5.448341 -3.849168 18.74917
------------------------------------------------------------------------------
diff = mean(AB) - mean(BA) t = 1.3674
H0: diff = 0 Degrees of freedom = 22
Ha: diff < 0 Ha: diff != 0 Ha: diff > 0
Pr(T < t) = 0.9073 Pr(|T| > |t|) = 0.1853 Pr(T > t) = 0.0927
.
. * design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
. * power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
. * (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
. scalar sig_s = 12
. scalar sig_e = 6
. scalar lam = -2
. count if grp == 1
12
. scalar n_ab = r(N)
. count if grp == 2
12
. scalar n_ba = r(N)
. scalar des_se = sqrt((4 * sig_s^2 + 2 * sig_e^2) * (1 / n_ab + 1 / n_ba))
. scalar des_df = n_ab + n_ba - 2
. scalar des_t = invttail(des_df, 0.025)
. scalar des_power = nt(des_df, lam / des_se, -des_t) + nttail(des_df, lam / des_se, des_t)
. scalar des_mde80 = (des_t + invttail(des_df, 0.20)) * des_se
. display "design SE of the carry-over estimate: " %7.4f des_se
design SE of the carry-over estimate: 10.3923
. display "power of the carry-over test at lambda = -2: " %6.4f des_power
power of the carry-over test at lambda = -2: 0.0539
. display "carry-over detected with 80% power: " %7.4f des_mde80
carry-over detected with 80% power: 30.4717R: การวิเคราะห์เดียวกันด้วย lme4 และ lmerTest
# Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any real drug)
# Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 12 participants, two periods.
# Packages: lme4 and lmerTest, plus pbkrtest for the Kenward-Roger degrees of freedom.
suppressPackageStartupMessages(library(lmerTest))
# read the simulated dataset (not published with this article)
dat <- read.csv("crossover.csv")
str(dat)
print(head(dat, 6))
# participants per sequence (one row per participant)
ids <- unique(dat[, c("id", "seq")])
print(table(ids$seq))
# cell means by sequence and period
print(round(tapply(dat$sbp, list(seq = dat$seq, period = dat$period), mean), 2))
# mixed model: treatment and period as fixed effects, a random intercept per participant (REML, Kenward-Roger df)
dat$trt <- factor(dat$trt)
dat$period <- factor(dat$period)
dat$carry <- factor(dat$carry)
mm <- lmer(sbp ~ trt + period + (1 | id), data = dat, REML = TRUE)
smm <- summary(mm, ddf = "Kenward-Roger")
print(smm)
cm <- coef(smm)
tq <- qt(0.975, cm["trt1", "df"])
# 95% CI of the treatment coefficient (A minus B)
print(round(cm["trt1", "Estimate"] + c(lo = -1, hi = 1) * tq * cm["trt1", "Std. Error"], 4))
# mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
mmc <- lmer(sbp ~ trt + period + carry + (1 | id), data = dat, REML = TRUE)
smmc <- summary(mmc, ddf = "Kenward-Roger")
print(smmc)
cc <- coef(smmc)
tqc <- qt(0.975, cc["trt1", "df"])
# 95% CI of the treatment coefficient in this model
print(round(cc["trt1", "Estimate"] + c(lo = -1, hi = 1) * tqc * cc["trt1", "Std. Error"], 4))
# one row per participant: period difference d = period 1 minus period 2, and the participant total
w <- reshape(dat[, c("id", "seq", "period", "sbp")], idvar = c("id", "seq"),
timevar = "period", direction = "wide")
w$d <- w$sbp.1 - w$sbp.2
w$tot <- w$sbp.1 + w$sbp.2
w$seq <- factor(w$seq, levels = c("AB", "BA"))
# treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
td <- t.test(d ~ seq, data = w, var.equal = TRUE)
print(td)
dbar_ab <- unname(td$estimate[1])
dbar_ba <- unname(td$estimate[2])
delta <- (dbar_ab - dbar_ba) / 2
delta_se <- td$stderr / 2
tc <- qt(0.975, unname(td$parameter))
cat(sprintf("treatment (A minus B): %.4f, SE %.4f, 95%% CI %.4f to %.4f\n",
delta, delta_se, delta - tc * delta_se, delta + tc * delta_se))
# period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
per <- -(dbar_ab + dbar_ba) / 2
cat(sprintf("period (2 minus 1): %.4f, 95%% CI %.4f to %.4f\n", per, per - tc * delta_se, per + tc * delta_se))
# carry-over test: compare participant totals between sequences (between-participant, so low power)
tt <- t.test(tot ~ seq, data = w, var.equal = TRUE)
print(tt)
# design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
# power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
# (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
sigma_s <- 12; sigma_e <- 6; lambda <- -2
n_ab <- sum(ids$seq == "AB"); n_ba <- sum(ids$seq == "BA")
des_se <- sqrt((4 * sigma_s^2 + 2 * sigma_e^2) * (1 / n_ab + 1 / n_ba))
des_df <- n_ab + n_ba - 2
des_t <- qt(0.975, des_df)
des_power <- pt(-des_t, des_df, lambda / des_se) + 1 - pt(des_t, des_df, lambda / des_se)
des_mde80 <- (des_t + qt(0.80, des_df)) * des_se
cat(sprintf("design SE of the carry-over estimate: %.4f\n", des_se))
cat(sprintf("power of the carry-over test at lambda = -2: %.4f\n", des_power))
cat(sprintf("carry-over detected with 80%% power: %.4f\n", des_mde80))> suppressPackageStartupMessages(library(lmerTest))
> dat <- read.csv("crossover.csv")
> str(dat)
'data.frame': 48 obs. of 7 variables:
$ id : int 1 1 2 2 3 3 4 4 5 5 ...
$ seq : chr "BA" "BA" "AB" "AB" ...
$ period : int 1 2 1 2 1 2 1 2 1 2 ...
$ treatment: chr "B" "A" "A" "B" ...
$ trt : int 0 1 1 0 1 0 1 0 1 0 ...
$ carry : int 0 0 0 1 0 1 0 1 0 1 ...
$ sbp : num 146 134 152 150 161 ...
> print(head(dat, 6))
id seq period treatment trt carry sbp
1 1 BA 1 B 0 0 146.2
2 1 BA 2 A 1 0 134.5
3 2 AB 1 A 1 0 152.4
4 2 AB 2 B 0 1 149.6
5 3 AB 1 A 1 0 160.7
6 3 AB 2 B 0 1 156.8
> ids <- unique(dat[, c("id", "seq")])
> print(table(ids$seq))
AB BA
12 12
> print(round(tapply(dat$sbp, list(seq = dat$seq, period = dat$period),
+ mean), 2))
period
seq 1 2
AB 145.72 149.60
BA 150.01 137.86
> dat$trt <- factor(dat$trt)
> dat$period <- factor(dat$period)
> dat$carry <- factor(dat$carry)
> mm <- lmer(sbp ~ trt + period + (1 | id), data = dat,
+ REML = TRUE)
> smm <- summary(mm, ddf = "Kenward-Roger")
> print(smm)
Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: sbp ~ trt + period + (1 | id)
Data: dat
REML criterion at convergence: 321
Scaled residuals:
Min 1Q Median 3Q Max
-2.0487 -0.4621 -0.1869 0.6003 1.4972
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 26.90 5.187
Residual 38.61 6.214
Number of obs: 48, groups: id, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 151.871 1.880 44.797 80.784 < 2e-16 ***
trt1 -8.017 1.794 22.000 -4.469 0.000192 ***
period2 -4.133 1.794 22.000 -2.304 0.031028 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trt1
trt1 -0.477
period2 -0.477 0.000
> cm <- coef(smm)
> tq <- qt(0.975, cm["trt1", "df"])
> print(round(cm["trt1", "Estimate"] + c(lo = -1, hi = 1) *
+ tq * cm["trt1", "Std. Error"], 4))
lo hi
-11.7367 -4.2966
> mmc <- lmer(sbp ~ trt + period + carry + (1 | id),
+ data = dat, REML = TRUE)
> smmc <- summary(mmc, ddf = "Kenward-Roger")
> print(smmc)
Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: sbp ~ trt + period + carry + (1 | id)
Data: dat
REML criterion at convergence: 313.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.1949 -0.5026 -0.0767 0.4785 1.5487
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 25.22 5.022
Residual 38.61 6.214
Number of obs: 48, groups: id, 24
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 150.008 2.306 38.059 65.041 <2e-16 ***
trt1 -4.292 3.262 38.059 -1.316 0.1961
period2 -7.858 3.262 38.059 -2.409 0.0209 *
carry1 7.450 5.448 22.000 1.367 0.1853
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trt1 perid2
trt1 -0.707
period2 0.279 -0.698
carry1 -0.591 0.835 -0.835
> cc <- coef(smmc)
> tqc <- qt(0.975, cc["trt1", "df"])
> print(round(cc["trt1", "Estimate"] + c(lo = -1, hi = 1) *
+ tqc * cc["trt1", "Std. Error"], 4))
lo hi
-10.8943 2.3110
> w <- reshape(dat[, c("id", "seq", "period", "sbp")],
+ idvar = c("id", "seq"), timevar = "period", direction = "wide")
> w$d <- w$sbp.1 - w$sbp.2
> w$tot <- w$sbp.1 + w$sbp.2
> w$seq <- factor(w$seq, levels = c("AB", "BA"))
> td <- t.test(d ~ seq, data = w, var.equal = TRUE)
> print(td)
Two Sample t-test
data: d by seq
t = -4.4691, df = 22, p-value = 0.0001918
alternative hypothesis: true difference in means between group AB and group BA is not equal to 0
95 percent confidence interval:
-23.473498 -8.593169
sample estimates:
mean in group AB mean in group BA
-3.883333 12.150000
> dbar_ab <- unname(td$estimate[1])
> dbar_ba <- unname(td$estimate[2])
> delta <- (dbar_ab - dbar_ba)/2
> delta_se <- td$stderr/2
> tc <- qt(0.975, unname(td$parameter))
> cat(sprintf("treatment (A minus B): %.4f, SE %.4f, 95%% CI %.4f to %.4f\n",
+ delta, delta_se, delta - tc * delta_se, delta + tc * delta_se))
treatment (A minus B): -8.0167, SE 1.7938, 95% CI -11.7367 to -4.2966
> per <- -(dbar_ab + dbar_ba)/2
> cat(sprintf("period (2 minus 1): %.4f, 95%% CI %.4f to %.4f\n",
+ per, per - tc * delta_se, per + tc * delta_se))
period (2 minus 1): -4.1333, 95% CI -7.8534 to -0.4133
> tt <- t.test(tot ~ seq, data = w, var.equal = TRUE)
> print(tt)
Two Sample t-test
data: tot by seq
t = 1.3674, df = 22, p-value = 0.1853
alternative hypothesis: true difference in means between group AB and group BA is not equal to 0
95 percent confidence interval:
-3.849168 18.749168
sample estimates:
mean in group AB mean in group BA
295.3167 287.8667
> sigma_s <- 12
> sigma_e <- 6
> lambda <- -2
> n_ab <- sum(ids$seq == "AB")
> n_ba <- sum(ids$seq == "BA")
> des_se <- sqrt((4 * sigma_s^2 + 2 * sigma_e^2) * (1/n_ab +
+ 1/n_ba))
> des_df <- n_ab + n_ba - 2
> des_t <- qt(0.975, des_df)
> des_power <- pt(-des_t, des_df, lambda/des_se) + 1 -
+ pt(des_t, des_df, lambda/des_se)
> des_mde80 <- (des_t + qt(0.8, des_df)) * des_se
> cat(sprintf("design SE of the carry-over estimate: %.4f\n",
+ des_se))
design SE of the carry-over estimate: 10.3923
> cat(sprintf("power of the carry-over test at lambda = -2: %.4f\n",
+ des_power))
power of the carry-over test at lambda = -2: 0.0539
> cat(sprintf("carry-over detected with 80%% power: %.4f\n",
+ des_mde80))
carry-over detected with 80% power: 30.4717ความเข้าใจผิดที่พบบ่อยและวิธีแก้
-
ยกการเปลี่ยนแปลงในช่วงเวลาที่ 2 ให้เป็นผลของยาในช่วงเวลาที่ 2
การเปลี่ยนแปลงจากช่วงเวลาที่ 1 ไปช่วงเวลาที่ 2 ภายในลำดับเดียวถูกอ่านเป็นผลของยาที่ได้ในช่วงเวลาที่ 2 ทั้งที่ส่วนหนึ่งอาจเป็นของช่วงเวลาที่ 2 เอง
วิธีแก้: เปรียบเทียบระหว่างลำดับ ไม่ใช่ระหว่างช่วงเวลาภายในลำดับเดียว ครึ่งหนึ่งของผลต่างระหว่างผลต่างระหว่างช่วงเวลาเฉลี่ยทั้งสองตัดการเลื่อนใด ๆ ที่ช่วงเวลาที่ 2 ได้รับร่วมกันออกไป
-
มองผลของลำดับเป็นผลแยกที่ประมาณได้
ในการศึกษาแบบไขว้ 2x2 ความต่างระหว่างกลุ่ม AB กับ BA แยกเป็นผลของลำดับ ผลตกค้าง และปฏิกิริยาสัมพันธ์ระหว่างการรักษากับช่วงเวลาไม่ได้
วิธีแก้: รายงานการเปรียบเทียบผลรวมเป็นปริมาณเดียวที่ปนกัน และป้องกันผลตกค้างด้วยการออกแบบ
-
"การทดสอบผลตกค้างที่ไม่มีนัยสำคัญแสดงว่าไม่มีผลตกค้าง"
ในการศึกษาจำลอง การทดสอบให้ P = 0.19 ทั้งที่สร้างผลตกค้าง -2 mmHg ไว้ในข้อมูล
วิธีแก้: การทดสอบผลตกค้างเปรียบเทียบผลรวมของผู้ร่วมวิจัยแต่ละคนระหว่างลำดับการได้รับยา ซึ่งเป็นการเปรียบเทียบระหว่างบุคคลที่มี power ต่ำ ผลที่ไม่มีนัยสำคัญจึงไม่ได้แสดงว่าไม่มีผลตกค้าง ผลตกค้างป้องกันได้ด้วยช่วงล้างยาที่เพียงพอ ไม่ใช่กำจัดได้ด้วยการทดสอบเบื้องต้น
-
ทดสอบผลตกค้างก่อน แล้วถอยไปใช้ช่วงเวลาที่ 1
ขั้นตอนสองขั้นทำให้ความคลาดเคลื่อนประเภทที่ 1 สูงเกินจริงและทำให้ค่าประมาณผลของการรักษาเอนเอียง [5]
วิธีแก้: กำหนดการวิเคราะห์ด้วยผลต่างระหว่างช่วงเวลา หรือแบบจำลองผสมที่เทียบเท่า ไว้ล่วงหน้า และวางแผนช่วงล้างยาที่เพียงพอ
สิ่งที่ควรทำในการวิเคราะห์ของคุณเอง
- วางแผนกำจัดผลตกค้างตั้งแต่ขั้นออกแบบ ด้วยภาวะที่คงที่ และช่วงล้างยาที่ยาวพอตามระยะเวลาออกฤทธิ์ที่ทราบของยา
- พิจารณากำหนดการวิเคราะห์หลักไว้ล่วงหน้าเป็นค่าประมาณจากผลต่างระหว่างช่วงเวลา หรือเมื่อข้อมูลครบถ้วนก็คือแบบจำลองผสมที่มีการรักษา ช่วงเวลา และจุดตัดแกนสุ่มรายผู้เข้าร่วม ซึ่งเทียบเท่ากัน
- รายงานค่าประมาณผลของช่วงเวลาไว้ข้างค่าประมาณผลของการรักษา เพื่อให้ผู้อ่านเห็นว่าเวลาเพียงอย่างเดียวขยับผลลัพธ์ไปเท่าไร
- หากรายงานค่าประมาณผลตกค้าง ให้นำเสนอว่าเป็นการเปรียบเทียบระหว่างผู้เข้าร่วมที่มีอำนาจการทดสอบต่ำ ซึ่งไม่ใช้กำหนดวิธีวิเคราะห์
- ตรวจรายงานกับส่วนขยายของ CONSORT (CONSORT extension) สำหรับการทดลองแบบไขว้ที่มีการสุ่ม ซึ่งเป็นแนวทางการรายงาน CONSORT ที่ปรับให้เหมาะกับแบบแผนนี้ [6] รายการหนึ่งในนั้นขอให้ระบุเหตุผลของการใช้แบบแผนไขว้ และวิธีทางสถิติที่คำนึงถึงลักษณะข้อมูลที่จับคู่กัน
อภิธานศัพท์
- crossover trial (การศึกษาแบบไขว้)
- การศึกษาที่ผู้เข้าร่วมแต่ละคนได้รับทุกการรักษาที่ศึกษา ช่วงเวลาละหนึ่งแบบ ตามลำดับที่สุ่มไว้
- sequence (ลำดับ)
- ลำดับที่ผู้เข้าร่วมได้รับการรักษา ในการศึกษาแบบไขว้ 2x2 คือ AB หรือ BA
- period effect (ผลของช่วงเวลา)
- การเลื่อนของผลลัพธ์ที่ผู้เข้าร่วมทุกคนได้รับร่วมกันในช่วงเวลาหนึ่งเมื่อเทียบกับอีกช่วงเวลา ไม่ว่าจะได้รับการรักษาแบบใด
- carry-over effect (ผลตกค้าง)
- ผลของการรักษาในช่วงเวลาก่อนหน้าที่ยังคงอยู่ในช่วงเวลาถัดมา
- washout period (ช่วงล้างยา)
- ช่วงที่ไม่มีการรักษาคั่นระหว่างช่วงเวลา เพื่อให้ผลของการรักษาแรกจางหายไป
- aliasing (การปนกันของผล)
- สองผลถือว่าปนกันเมื่อแบบแผนทำให้ทั้งสองมีรูปแบบเดียวกันในข้อมูล จึงไม่มีการวิเคราะห์ใดแยกได้
- within-participant comparison (การเปรียบเทียบภายในผู้เข้าร่วม)
- การเปรียบเทียบค่าที่วัดจากคนเดียวกัน ความต่างที่คงที่ระหว่างคนจึงหักล้างกัน
- power (อำนาจการทดสอบ)
- ความน่าจะเป็นที่การทดสอบจะประกาศว่ามีผล เมื่อมีผลขนาดที่ระบุอยู่จริง
- CONSORT extension (ส่วนขยายของ CONSORT)
- แนวทางการรายงาน CONSORT สำหรับการทดลองแบบสุ่มฉบับที่ปรับให้เหมาะกับแบบแผนหนึ่ง ฉบับสำหรับการศึกษาแบบไขว้ตีพิมพ์ในปี 2019
เอกสารอ้างอิง
- Senn S. Cross-over trials in clinical research. 2nd ed. Wiley; 2002. https://doi.org/10.1002/0470854596
- Jones B, Kenward MG. Design and analysis of cross-over trials. 3rd ed. Chapman and Hall/CRC; 2014. https://doi.org/10.1201/b17537
- Hills M, Armitage P. The two-period cross-over clinical trial. Br J Clin Pharmacol. 1979;8(1):7-20. https://doi.org/10.1111/j.1365-2125.1979.tb05903.x
- Grizzle JE. The two-period change-over design and its use in clinical trials. Biometrics. 1965;21(2):467-80. https://doi.org/10.2307/2528104
- Freeman PR. The performance of the two-stage analysis of two-treatment, two-period crossover trials. Stat Med. 1989;8(12):1421-32. https://doi.org/10.1002/sim.4780081202
- Dwan K, Li T, Altman DG, Elbourne D. CONSORT 2010 statement: extension to randomised crossover trials. BMJ. 2019;366:l4378. https://doi.org/10.1136/bmj.l4378
ประเด็นสำคัญ
- ในการศึกษาแบบไขว้ 2x2 ครึ่งหนึ่งของความต่างระหว่างค่าเฉลี่ยของผลต่างระหว่างช่วงเวลาของสองลำดับประมาณผลของการรักษา และผลของช่วงเวลาใด ๆ ถูกตัดออก
- เมื่อข้อมูลครบถ้วน แบบจำลองผสมที่มีจุดตัดแกนสุ่ม พร้อมการรักษาและช่วงเวลา ประมาณครึ่งหนึ่งของความต่างระหว่างค่าเฉลี่ยของผลต่างระหว่างช่วงเวลา โดยได้ค่า P และองศาอิสระเดียวกับการทดสอบทีแบบความแปรปรวนรวม
- ผลตกค้างปนอยู่กับลำดับ มันเลื่อนค่าประมาณผลของการรักษาไปลบครึ่งหนึ่งของขนาดของมัน และประมาณได้เฉพาะจากการเปรียบเทียบระหว่างผู้เข้าร่วม
- การทดสอบผลตกค้างมีอำนาจการทดสอบต่ำ ผลที่ไม่มีนัยสำคัญจึงไม่ได้แสดงว่าไม่มีผลตกค้าง
- ป้องกันผลตกค้างด้วยช่วงล้างยาที่เพียงพอ และอย่าให้การทดสอบผลตกค้างเบื้องต้นเป็นตัวกำหนดวิธีวิเคราะห์
อ่านต่อในวิกิ: [[treatment-by-time-interaction-longitudinal-trials-th]] [[crossover-trials-design-analysis-guide]] [[mixed-model-series-guide-th]]