Causal Mediation: ทำไม a x b จึงใช้ไม่ได้อีกต่อไป

Read the English version

บทคัดย่อ

การวิเคราะห์ตัวแปรสื่อกลาง (mediation analysis) ถามว่าผลของสิ่งสัมผัส (exposure) ต่อผลลัพธ์ผ่านตัวแปรสื่อกลาง (mediator) ซึ่งเป็นตัวแปรบนเส้นทางเหตุและผลระหว่างทั้งสองมากเพียงใด สูตรที่คุ้นเคยคือการคูณความชันสองค่า ผลคูณ a x b เท่ากับ natural indirect effect ก็ต่อเมื่อแบบจำลองทั้งสองเป็นเชิงเส้นและไม่มีปฏิกิริยาสัมพันธ์ระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง และจะมีความหมายเชิงเหตุผลได้ก็ต่อเมื่อไม่มีตัวกวนที่ไม่ได้วัดระหว่างสิ่งสัมผัสกับผลลัพธ์ ระหว่างตัวแปรสื่อกลางกับผลลัพธ์ และระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง รวมทั้งไม่มีตัวกวนระหว่างตัวแปรสื่อกลางกับผลลัพธ์ที่ถูกสิ่งสัมผัสส่งผลกระทบ ใน cohort จำลองของผู้ใหญ่ 2,000 คน การออกกำลังกายลดความดันโลหิตส่วนหนึ่งผ่านการเปลี่ยนแปลงของน้ำหนัก และน้ำหนักแต่ละกิโลกรัมมีผลมากขึ้นเมื่อออกกำลังกายมากขึ้น สำหรับการออกกำลังกาย 2.5 ชั่วโมงต่อสัปดาห์เทียบกับ 0 ชั่วโมง ผลทางอ้อมตามธรรมชาติที่แท้จริงในประชากรจำลองคือ -4.95 mmHg ขณะที่วิธีผลคูณมีเป้าหมายอยู่ที่ -4.06 mmHg บทความนี้นิยามผลเหล่านี้ด้วยผลลัพธ์ที่อาจเกิดขึ้น (potential outcomes) ให้สูตรในรูปปิดที่คงปฏิกิริยาสัมพันธ์ไว้ แสดงโค้ด Stata และ R พร้อมช่วงเชื่อมั่นจากวิธีเดลตา (delta method) ครอบคลุมผลลัพธ์แบบสองค่าที่พบบ่อย และระบุข้อสมมติข้ามโลก (cross-world assumption)

🔍 คลิกเพื่อขยาย คลิกที่ใดก็ได้เพื่อปิด

ผลคูณของความชันสองค่า และคำถามของผู้ทบทวนร่างรายงาน

ทีมสุขภาพชุมชนทีมหนึ่งจัดโปรแกรมออกกำลังกายมาเป็นเวลาหนึ่งปี ในผู้ใหญ่ 2,000 คน ผู้ที่ออกกำลังกายมากกว่ามีความดันโลหิตซิสโตลิก (systolic blood pressure, SBP) ต่ำกว่าเมื่อสิ้นปี และน้ำหนักลดลงมากกว่าด้วย ทีมต้องการทราบว่าความดันที่ลดลงเกิดผ่านการลดน้ำหนักเป็นสัดส่วนเท่าใด โปรแกรมนี้เป็นของสมมติและข้อมูลเป็นข้อมูลจำลอง

นักวิเคราะห์ใช้สูตรที่คุ้นเคย คือถดถอยการเปลี่ยนแปลงของน้ำหนักบนการออกกำลังกาย ถดถอยความดันโลหิตบนการออกกำลังกายและการเปลี่ยนแปลงของน้ำหนัก แล้วคูณความชันสองค่านั้น ผลคูณถูกรายงานเป็นผลทางอ้อม (indirect effect) และความชันของการออกกำลังกายในแบบจำลองที่สองถูกรายงานเป็นผลทางตรง (direct effect)

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

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

วิธีผลคูณและกรณีที่ใช้ได้

ให้ $A$ แทนสิ่งสัมผัส (exposure) คือชั่วโมงออกกำลังกายต่อสัปดาห์ และ $M$ แทนตัวแปรสื่อกลาง คือการเปลี่ยนแปลงของน้ำหนักเป็นกิโลกรัมตลอดหนึ่งปี โดยค่าติดลบหมายถึงน้ำหนักลด $Y$ คือผลลัพธ์ ได้แก่ SBP หน่วย mmHg และ $C$ แทนตัวแปรร่วม ณ จุดเริ่มต้น ในที่นี้คืออายุและเพศ ซึ่งอาจเป็นตัวกวนของความสัมพันธ์เหล่านี้ ตัวแปรสื่อกลาง (mediator) คือตัวแปรบนเส้นทางเหตุและผลจากสิ่งสัมผัสไปสู่ผลลัพธ์ การออกกำลังกายเปลี่ยนน้ำหนัก และน้ำหนักเปลี่ยนความดัน

วิธีถดถอยแบบดั้งเดิมในการวิเคราะห์ตัวแปรสื่อกลางประมาณแบบจำลองเชิงเส้นสองแบบ แบบหนึ่งสำหรับตัวแปรสื่อกลางและอีกแบบสำหรับผลลัพธ์ [1] แบบจำลองของตัวแปรสื่อกลางถดถอย $M$ บน $A$ และ $C$ ความชันของการออกกำลังกายในแบบจำลองนี้คือ $\beta_1$ ซึ่งเป็นการเปลี่ยนแปลงของค่าเฉลี่ยการเปลี่ยนแปลงน้ำหนักต่อชั่วโมงต่อสัปดาห์ที่เพิ่มขึ้น แบบจำลองของผลลัพธ์ถดถอย $Y$ บน $A$, $M$ และ $C$ โดยไม่มีพจน์ผลคูณ ความชันของน้ำหนักในแบบจำลองนี้คือ $b$ ซึ่งเป็นการเปลี่ยนแปลงของค่าเฉลี่ย SBP ต่อกิโลกรัม และความชันของการออกกำลังกายคือ $c'$ ซึ่งอ่านเป็นผลทางตรงต่อชั่วโมงต่อสัปดาห์

ในสูตรนี้ $a^*$ และ $a$ คือระดับการออกกำลังกายสองระดับที่นำมาเปรียบเทียบ โดย $a^* = 0$ และ $a = 2.5$ ชั่วโมงต่อสัปดาห์ ความกว้างของคอนทราสต์ (contrast) $a - a^*$ จึงเท่ากับ 2.5 ชั่วโมง ผลทางตรงในคอนทราสต์เดียวกันคือ $c' \times (a - a^*)$ ในสูตร a x b ที่คุ้นเคย ความชันจากสิ่งสัมผัสไปยังตัวแปรสื่อกลางถูกเรียกว่า $a$ แต่บทความนี้สงวน $a$ ไว้เป็นระดับของสิ่งสัมผัส และเขียนความชันนั้นเป็น $\beta_1$

ผลคูณ a x b เท่ากับ natural indirect effect ก็ต่อเมื่อแบบจำลองทั้งสองเป็นเชิงเส้นและไม่มีปฏิกิริยาสัมพันธ์ระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง และจะมีความหมายเชิงเหตุผลได้ก็ต่อเมื่อไม่มีตัวกวนที่ไม่ได้วัดระหว่างสิ่งสัมผัสกับผลลัพธ์ ระหว่างตัวแปรสื่อกลางกับผลลัพธ์ และระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง รวมทั้งไม่มีตัวกวนระหว่างตัวแปรสื่อกลางกับผลลัพธ์ที่ถูกสิ่งสัมผัสส่งผลกระทบ [2] natural indirect effect (ผลทางอ้อมตามธรรมชาติ) มีนิยามอย่างเป็นทางการในหัวข้อถัดไป ในกรณีเชิงเส้นที่ไม่มีปฏิกิริยาสัมพันธ์บนมาตรวัดผลต่างหน่วย mmHg การออกกำลังกายเลื่อนค่าเฉลี่ยของน้ำหนักไป $\beta_1$ ต่อชั่วโมงต่อสัปดาห์ น้ำหนักทุกกิโลกรัมเลื่อนค่าเฉลี่ยของความดันไป $b$ เท่ากันหมด และการเลื่อนทั้งสองคูณกันได้ เครื่องหมายสำคัญพอ ๆ กับขนาด ดังที่ตัวอย่างคำนวณด้วยมือแสดง

วิธีผลคูณใน cohort จำลอง

cohort จำลองที่ใช้ตลอดบทความนี้มีผู้ใหญ่ 2,000 คน อายุ 30 ถึง 80 ปี แต่ละคนมีข้อมูลการออกกำลังกายรายสัปดาห์ การเปลี่ยนแปลงของน้ำหนักตลอดหนึ่งปี SBP เมื่อครบหนึ่งปี อายุ และเพศ เครื่องหมายหมวก (hat) แสดงค่าประมาณ เมื่อประมาณแบบจำลองเชิงเส้นสองแบบโดยใช้อายุและเพศเป็นตัวแปรร่วม จะได้ $\hat\beta_1 = -1.19$ กิโลกรัมต่อชั่วโมงต่อสัปดาห์ และ $\hat b = 1.37$ mmHg ต่อกิโลกรัม ช่วงเชื่อมั่น 95% ที่พิมพ์ในผลลัพธ์สองชุดด้านล่างครอบคลุมค่าที่ความชันเหล่านี้มีในประชากรจำลอง

ในคอนทราสต์ 2.5 ชั่วโมง ผลคูณของความชันสองค่ากับความกว้างของคอนทราสต์คือ -4.09 mmHg ช่วงเชื่อมั่น 95% จากวิธีเดลตา (delta method) คือ -4.81 ถึง -3.37 mmHg วิธีเดลตาให้ค่าคลาดเคลื่อนมาตรฐานของฟังก์ชันของสัมประสิทธิ์ที่ประมาณได้ เช่นผลคูณ โดยถือว่าฟังก์ชันนั้นเกือบเป็นเส้นตรงในบริเวณใกล้ค่าประมาณ สำหรับผลคูณของความชันสองค่า วิธีนี้ยังเรียกว่าค่าคลาดเคลื่อนมาตรฐานของ Sobel ทุกช่วงเชื่อมั่นในตารางผลลัพธ์สองตารางเป็นช่วงเชื่อมั่น 95% จากวิธีเดลตา ส่วนผลลัพธ์ของคำสั่ง mediate ของ Stata และแพ็กเกจ mediation ของ R ใช้วิธีหาช่วงอื่น ตามที่คำอธิบายใต้ผลลัพธ์ระบุไว้

Stata: แบบจำลองของตัวแปรสื่อกลาง

* Mediator model: weight change on exercise and the confounders (age, sex)
regress wtchg exercise age male

. * Mediator model: weight change on exercise and the confounders (age, sex)
. regress wtchg exercise age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(3, 1996)      =    412.53
       Model |  7937.34882         3  2645.78294   Prob > F        =    0.0000
    Residual |  12801.4031     1,996  6.41352861   R-squared       =    0.3827
-------------+----------------------------------   Adj R-squared   =    0.3818
       Total |  20738.7519     1,999  10.3745632   Root MSE        =    2.5325

------------------------------------------------------------------------------
       wtchg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.194437   .0435007   -27.46   0.000    -1.279749   -1.109126
         age |   .0236362   .0044603     5.30   0.000     .0148888    .0323836
        male |   .2156961   .1139233     1.89   0.058     -.007725    .4391173
       _cons |  -.1197085   .3035733    -0.39   0.693    -.7150622    .4756452
------------------------------------------------------------------------------

Stata: แบบจำลองของผลลัพธ์ที่ไม่มีพจน์ปฏิกิริยาสัมพันธ์

* Product method: outcome model without the interaction gives b, then a x b over the 2.5-hour contrast
regress sbp exercise wtchg age male

. * Product method: outcome model without the interaction gives b, then a x b over the 2.5-hour contrast
. regress sbp exercise wtchg age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(4, 1995)      =    403.71
       Model |  263039.698         4  65759.9245   Prob > F        =    0.0000
    Residual |  324960.139     1,995  162.887288   R-squared       =    0.4473
-------------+----------------------------------   Adj R-squared   =    0.4462
       Total |  587999.837     1,999  294.146992   Root MSE        =    12.763

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.979195   .2573194    -7.69   0.000    -2.483838   -1.474552
       wtchg |   1.369749   .1128015    12.14   0.000     1.148528     1.59097
         age |   .4485983   .0226358    19.82   0.000      .404206    .4929905
        male |   3.890667   .5746421     6.77   0.000     2.763706    5.017629
       _cons |   103.2818   1.529944    67.51   0.000     100.2814    106.2823
------------------------------------------------------------------------------

สามกรณีที่วิธีผลคูณใช้ไม่ได้

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

ในสมการนี้ $E[Y \mid a, m, c]$ คือ SBP เฉลี่ยที่ระดับการออกกำลังกาย $a$ การเปลี่ยนแปลงของน้ำหนัก $m$ และตัวแปรร่วม $c$ สัมประสิทธิ์ $\theta_1$ คือความชันของการออกกำลังกายเมื่อน้ำหนักไม่เปลี่ยน $\theta_2$ คือความชันของน้ำหนักเมื่อไม่ออกกำลังกาย และ $\theta_4$ คือสัมประสิทธิ์ของตัวแปรร่วม ส่วน $\theta_3$ คือปฏิกิริยาสัมพันธ์ระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง (exposure-mediator interaction) หมายถึงความชันของน้ำหนักเปลี่ยนไปเท่าใดต่อชั่วโมงออกกำลังกายต่อสัปดาห์ที่เพิ่มขึ้น ความชันของน้ำหนักที่ระดับการออกกำลังกาย $a$ จึงเท่ากับ $\theta_2 + \theta_3 a$

cohort จำลองสร้างขึ้นโดยใช้ $\theta_2 = 0.4$ และ $\theta_3 = 0.5$ น้ำหนักที่เปลี่ยนไปทุกกิโลกรัมจึงเปลี่ยน SBP ไป 0.4 mmHg เมื่อไม่ออกกำลังกาย และเปลี่ยนไป $0.4 + 0.5 \times 2.5 = 1.65$ mmHg เมื่อออกกำลังกาย 2.5 ชั่วโมงต่อสัปดาห์ ตัวเลขทั้งสองไม่ใช่ $b$ ของวิธีผลคูณ แบบจำลองที่ไม่มีพจน์ผลคูณผสมความชันของน้ำหนักข้ามระดับการออกกำลังกายต่าง ๆ ใน cohort และในตัวอย่างจำลองขนาดใหญ่มาก $b$ จะลงตัวที่ 1.35 mmHg ต่อกิโลกรัม

การผสมนั้นคือความล้มเหลวประการแรก เป้าหมายของวิธีผลคูณในประชากรจำลองคือ -4.06 mmHg ขณะที่ natural indirect effect ซึ่งนิยามในหัวข้อถัดไปคือ -4.95 mmHg ช่องว่างเกิดจากการใช้ความชันของน้ำหนักที่เฉลี่ยแล้วเพียงค่าเดียว ทั้งที่ผลทางอ้อมต้องใช้ความชันที่ระดับการออกกำลังกายซึ่งนำมาเปรียบเทียบ

ความล้มเหลวประการที่สองคือแบบจำลองของผลลัพธ์ที่ไม่เป็นเชิงเส้น สำหรับผลลัพธ์แบบใช่หรือไม่ใช่ เช่นความดันโลหิตสูง แบบจำลองของผลลัพธ์มักเป็นการถดถอยโลจิสติก ซึ่งสัมประสิทธิ์คือการเปลี่ยนแปลงของ log odds ผลคูณของ $\beta_1$ กับสัมประสิทธิ์ log odds ไม่ใช่ผลบนมาตรวัดผลต่างความเสี่ยง (risk-difference scale) เมื่อคูณด้วยความกว้างของคอนทราสต์แล้วยกกำลังด้วยเลขชี้กำลัง ผลคูณนั้นประมาณ natural indirect effect บนมาตรวัด odds ratio ได้ แต่เฉพาะเมื่อผลลัพธ์พบน้อยและแบบจำลองโลจิสติกไม่มีพจน์ผลคูณระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง นั่นคือไม่มีปฏิกิริยาสัมพันธ์บนมาตรวัด log odds

ความล้มเหลวประการที่สามคือตัวแปรสื่อกลางแบบสองค่า เช่นการเข้าร่วมแผนลดเกลือระหว่างโปรแกรมหรือไม่ แบบจำลองของมันก็เป็นโลจิสติกเช่นกัน และ $\beta_1$ กลายเป็นการเปลี่ยนแปลงของ log odds ของตัวแปรสื่อกลาง ไม่ใช่การเลื่อนของค่าเฉลี่ย ผลทางอ้อมจึงขึ้นกับว่าความน่าจะเป็นของตัวแปรสื่อกลางขยับไปไกลเท่าใด ระยะนั้นขึ้นกับตัวแปรร่วมของแต่ละคน เพราะ log odds ที่เปลี่ยนเท่ากันขยับความน่าจะเป็นได้ไม่เท่ากันเมื่อความน่าจะเป็นตั้งต้นต่างกัน

เมื่อประมาณแบบจำลองของผลลัพธ์ที่มีพจน์ผลคูณกับ cohort จำลอง จะได้ $\hat\theta_3 = 0.55$ mmHg ต่อกิโลกรัมต่อชั่วโมงต่อสัปดาห์ โดย $\hat\theta_1 = -0.84$ mmHg ต่อชั่วโมงต่อสัปดาห์ และ $\hat\theta_2 = 0.29$ mmHg ต่อกิโลกรัม ช่วงเชื่อมั่น 95% ในผลลัพธ์ด้านล่างครอบคลุมค่าที่ใช้สร้างข้อมูล คือ 0.5, -1 และ 0.4 ค่าประมาณแต่ละค่าจึงอยู่ในขอบเขตความคลาดเคลื่อนจากการสุ่มของค่าที่ใช้สร้างมัน

Stata: แบบจำลองของผลลัพธ์ที่มีปฏิกิริยาสัมพันธ์ระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง

* Outcome model with the exposure-mediator interaction (theta1 A + theta2 M + theta3 A M)
regress sbp exercise wtchg exw age male

. * Outcome model with the exposure-mediator interaction (theta1 A + theta2 M + theta3 A M)
. regress sbp exercise wtchg exw age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(5, 1994)      =    358.74
       Model |  278452.465         5  55690.4931   Prob > F        =    0.0000
    Residual |  309547.372     1,994  155.239404   R-squared       =    0.4736
-------------+----------------------------------   Adj R-squared   =    0.4722
       Total |  587999.837     1,999  294.146992   Root MSE        =     12.46

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -.8351658   .2762007    -3.02   0.003    -1.376838   -.2934935
       wtchg |   .2898677   .1545066     1.88   0.061    -.0131437    .5928791
         exw |   .5457527   .0547717     9.96   0.000     .4383369    .6531686
         age |   .4668235   .0221736    21.05   0.000     .4233377    .5103092
        male |   3.880244   .5609905     6.92   0.000     2.780054    4.980433
       _cons |   101.6947   1.502064    67.70   0.000     98.74896    104.6405
------------------------------------------------------------------------------

ตั้งชื่อผลด้วยผลลัพธ์ที่อาจเกิดขึ้น

ผลลัพธ์ที่อาจเกิดขึ้น (potential outcome) คือค่าที่ผลลัพธ์จะมี หากกำหนดสิ่งสัมผัส และอาจรวมถึงตัวแปรสื่อกลาง ไว้ที่ระดับหนึ่ง [3, 4] ให้ $M^{a}$ แทนการเปลี่ยนแปลงของน้ำหนักที่บุคคลหนึ่งจะมี หากกำหนดการออกกำลังกายไว้ที่ $a$ และให้ $Y^{a,m}$ แทน SBP ที่บุคคลนั้นจะมี หากกำหนดการออกกำลังกายไว้ที่ $a$ และการเปลี่ยนแปลงของน้ำหนักไว้ที่ $m$ สำหรับบุคคลหนึ่งคนจะสังเกตเห็นค่าเหล่านี้ได้อย่างมากเพียงค่าเดียว นิยามด้านล่างจึงเปรียบเทียบค่าเฉลี่ย เขียนเป็น $E[\cdot]$

เมื่อซ้อนทั้งสองเข้าด้วยกันจะได้ $Y^{a,M^{a^*}}$ คือ SBP เมื่อการออกกำลังกายอยู่ที่ $a$ แต่การเปลี่ยนแปลงของน้ำหนักอยู่ที่ระดับที่จะเกิดขึ้นหากการออกกำลังกายอยู่ที่ $a^*$ เมื่อการออกกำลังกายอยู่ที่ $a$ การเปลี่ยนแปลงของน้ำหนักจะมีค่าตามธรรมชาติ $M^{a}$ ดังนั้น $Y^{a} = Y^{a,M^{a}}$ ผลตามธรรมชาติแต่ละตัวด้านล่างเปลี่ยนทีละอย่างเดียว

ผลทางตรงแบบควบคุม (controlled direct effect) กำหนดการเปลี่ยนแปลงของน้ำหนักไว้ที่ค่า $m$ เดียวกันสำหรับทุกคน และเปลี่ยนเฉพาะการออกกำลังกาย มันตอบว่าการออกกำลังกายที่เพิ่มขึ้นจะทำอะไรได้ หากสามารถตรึงการเปลี่ยนแปลงของน้ำหนักไว้ที่ $m$ ด้วยการแทรกแซง

ผลทางตรงตามธรรมชาติ (natural direct effect) เปลี่ยนการออกกำลังกายจาก $a^*$ เป็น $a$ แต่คงการเปลี่ยนแปลงของน้ำหนักของแต่ละคนไว้ที่ระดับที่จะเกิดขึ้นหากไม่ได้ออกกำลังกายเพิ่ม จึงเปิดเฉพาะเส้นทางที่อ้อมน้ำหนักไป

ผลทางอ้อมตามธรรมชาติ (natural indirect effect) คงการออกกำลังกายไว้ที่ $a$ และเลื่อนการเปลี่ยนแปลงของน้ำหนักของแต่ละคนจากระดับภายใต้ $a^*$ ไปเป็นระดับภายใต้ $a$ จึงเปิดเฉพาะเส้นทางที่ผ่านน้ำหนัก

ผลรวมทั้งหมด (total effect) เปรียบเทียบการออกกำลังกายสองระดับโดยให้น้ำหนักตอบสนองได้อย่างอิสระ ผลตามธรรมชาติทั้งสองบวกกันได้เท่ากับผลรวมทั้งหมด เพราะพจน์กลาง $Y^{a,M^{a^*}}$ ปรากฏในทั้งสองและหักล้างกัน

สัดส่วนของผลที่ผ่านตัวแปรสื่อกลาง (proportion mediated) คือ $\mathrm{NIE}/\mathrm{TE}$ ผลทางอ้อมในรูปสัดส่วนของผลรวมทั้งหมดบนมาตรวัดเดียวกัน หัวข้อ 'สัดส่วนของผลที่ผ่านตัวแปรสื่อกลาง ต้องอ่านด้วยความระมัดระวัง' แสดงว่าทำไมต้องระวัง ในแบบจำลองเชิงเส้นที่ไม่มีปฏิกิริยาสัมพันธ์บนมาตรวัด mmHg ผลทางตรงแบบควบคุมมีค่าเท่ากันทุกค่าของ $m$ และเท่ากับผลทางตรงตามธรรมชาติ เมื่อมีปฏิกิริยาสัมพันธ์ ทั้งสองมักต่างกัน

ปริมาณ $Y^{a,M^{a^*}}$ ผสมสองโลกเข้าด้วยกัน คือการออกกำลังกายตามโลกหนึ่งและการเปลี่ยนแปลงของน้ำหนักตามอีกโลกหนึ่ง จึงไม่เคยสังเกตเห็นในใครเลย นั่นคือสิ่งที่ทำให้ผลตามธรรมชาติบอกกลไกได้ และเป็นสิ่งเดียวกันที่ทำให้ระบุค่าได้ยาก (identify) คือคำนวณจากข้อมูลที่สังเกตได้ภายใต้ข้อสมมติที่ระบุไว้ (ดูหัวข้อ 'สิ่งที่ต้องมีจึงจะอ่านผลเหล่านี้เชิงเหตุผลได้')

สูตรในรูปปิดเมื่อการออกกำลังกายเปลี่ยนผลของน้ำหนัก

เมื่อแบบจำลองของตัวแปรสื่อกลางเป็นเชิงเส้น และแบบจำลองของผลลัพธ์เป็นเชิงเส้นที่คงพจน์ผลคูณไว้ ผลแต่ละตัวมีสูตรในรูปปิด (closed form) คือสูตรที่เขียนด้วยสัมประสิทธิ์ของแบบจำลอง [2, 5] แบบจำลองของตัวแปรสื่อกลางคือ

โดย $\beta_0$ คือจุดตัดแกน (intercept) $\beta_1$ คือความชันของการออกกำลังกาย และ $\beta_2$ คือสัมประสิทธิ์ของตัวแปรร่วม เมื่อใช้ร่วมกับแบบจำลองของผลลัพธ์ข้างต้น ผลสำหรับคอนทราสต์จาก $a^*$ ไป $a$ ที่ค่าตัวแปรร่วม $c$ เป็นดังนี้

ผลทางตรงแบบควบคุมขึ้นกับระดับ $m$ ที่ตรึงการเปลี่ยนแปลงของน้ำหนักไว้ เพราะความชันของการออกกำลังกายที่ระดับนั้นคือ $\theta_1 + \theta_3 m$

ผลทางตรงตามธรรมชาติคือผลทางตรงแบบควบคุมที่ประเมินที่ค่าเฉลี่ยของการเปลี่ยนแปลงน้ำหนักภายใต้ $a^*$ ซึ่งคือ $\beta_0 + \beta_1 a^* + \beta_2' c$

ผลทางอ้อมตามธรรมชาติคูณการเลื่อนของค่าเฉลี่ยการเปลี่ยนแปลงน้ำหนัก $\beta_1 (a - a^*)$ ด้วยความชันของน้ำหนักที่ระดับการออกกำลังกาย $a$ ซึ่งคือ $\theta_2 + \theta_3 a$

เมื่อกำหนด $\theta_3 = 0$ สูตรทั้งหมดจะยุบเหลือวิธีผลคูณ ผลทางตรงตามธรรมชาติกลายเป็น $\theta_1 (a - a^*)$ และผลทางอ้อมตามธรรมชาติกลายเป็น $\theta_2 \beta_1 (a - a^*)$ โดย $\theta_2$ เท่ากับ $b$ เมื่อ $\theta_3 \ne 0$ ผลทางอ้อมใช้ความชันของน้ำหนักที่ $a$ และผลทางตรงใช้ค่าเฉลี่ยการเปลี่ยนแปลงน้ำหนักที่ $a^*$ เนื่องจากแบบจำลองทั้งสองเป็นเชิงเส้น การแทนค่าเฉลี่ยของตัวแปรร่วมจึงให้ผลเฉลี่ยทั้ง cohort

cohort จำลอง: ค่าประมาณเทียบกับค่าจริงในข้อมูลจำลอง

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

ทุกช่วงเชื่อมั่นจากวิธีเดลตาในตารางครอบคลุมค่าจริงในข้อมูลจำลองของมัน ในตัวอย่างนี้ ผลทางอ้อมตามธรรมชาติอยู่ใกล้ค่าจริงในข้อมูลจำลอง คือ -4.94 mmHg เทียบกับ -4.95 ส่วนผลทางตรงตามธรรมชาติอยู่ไกลกว่า คือ -0.32 mmHg เทียบกับ -1.00 มันเป็นผลต่างเล็ก ๆ ของพจน์ที่ประมาณได้สองพจน์ซึ่งใหญ่กว่า ช่วงเชื่อมั่นจากวิธีเดลตาของมันคือ -1.87 ถึง 1.23 mmHg จึงกว้าง

วิธีผลคูณ เมื่อใช้กับข้อมูลชุดเดียวกัน ให้ -4.09 mmHg โดยช่วงเชื่อมั่นจากวิธีเดลตาคือ -4.81 ถึง -3.37 ช่วงนี้ครอบคลุมเป้าหมายของวิธีผลคูณเอง คือ -4.06 และไม่ครอบคลุมผลทางอ้อมตามธรรมชาติที่แท้จริงในข้อมูลจำลอง คือ -4.95 ตัวอย่างที่ใหญ่ขึ้นจะทำให้ช่วงแคบลงรอบเป้าหมายที่ผิดเท่านั้น

ส่วนที่แสดงโค้ดด้านล่างแสดงสูตรในรูปปิดเป็นโค้ด คำสั่ง nlcom ของ Stata ซึ่งย่อมาจากการรวมแบบไม่เป็นเชิงเส้นของค่าประมาณ (nonlinear combinations of estimates) ใช้วิธีเดลตาเมื่อซ้อนแบบจำลองสองแบบเข้าด้วยกันแล้ว ส่วนใน R ใช้อนุพันธ์เชิงตัวเลขของสูตรเดียวกัน คำสั่ง mediate ที่มีในตัว Stata คำสั่ง paramed ที่ผู้ใช้เขียนขึ้น และแพ็กเกจ mediation ของ R ซึ่งใช้วิธีจำลองของ Imai และคณะ [6] ให้ค่าประมาณเท่ากันหรือใกล้เคียงมาก

ผลลัพธ์ต่อเนื่อง: ค่าจริงในข้อมูลจำลอง ค่าประมาณ และช่วงเชื่อมั่นจากวิธีเดลตา

Stata: สูตรในรูปปิดและช่วงเชื่อมั่นจากวิธีเดลตาด้วย nlcom

* Stack both models (block-diagonal variance) so nlcom can apply the delta method
matrix coleq bm = med
matrix coleq by = out
matrix bb = bm, by
local names : colfullnames bb
matrix VV = (Vm, J(4, 6, 0) \ J(6, 4, 0), Vy)
matrix colnames VV = `names'
matrix rownames VV = `names'
ereturn post bb VV

* Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
local m0 "(_b[med:_cons] + _b[med:age]*(`agebar') + _b[med:male]*(`malebar'))"
local nde "((_b[out:exercise] + _b[out:exw]*`m0')*2.5)"
local nie "((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)"
local te "(`nde' + `nie')"
nlcom (nde: `nde') (nie: `nie') (te: `te') (pm: `nie'/`te') (cde0: _b[out:exercise]*2.5) (pnie: _b[out:wtchg]*_b[med:exercise]*2.5) (tnde: (_b[out:exercise] + _b[out:exw]*(`m0' + _b[med:exercise]*2.5))*2.5)

. * Stack both models (block-diagonal variance) so nlcom can apply the delta method
. matrix coleq bm = med

. matrix coleq by = out

. matrix bb = bm, by

. local names : colfullnames bb

. matrix VV = (Vm, J(4, 6, 0) \ J(6, 4, 0), Vy)

. matrix colnames VV = `names'

. matrix rownames VV = `names'

. ereturn post bb VV

.
. * Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
. local m0 "(_b[med:_cons] + _b[med:age]*(`agebar') + _b[med:male]*(`malebar'))"

. local nde "((_b[out:exercise] + _b[out:exw]*`m0')*2.5)"

. local nie "((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)"

. local te "(`nde' + `nie')"

. nlcom (nde: `nde') (nie: `nie') (te: `te') (pm: `nie'/`te') (cde0: _b[out:exercise]*2.5) (pnie: _b[out:wtchg]*_b[med:e
> xercise]*2.5) (tnde: (_b[out:exercise] + _b[out:exw]*(`m0' + _b[med:exercise]*2.5))*2.5)

         nde: ((_b[out:exercise] + _b[out:exw]*(_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5)
         nie: ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)
          te: (((_b[out:exercise] + _b[out:exw]*(_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5) + ((_
> b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5))
          pm: ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)/(((_b[out:exercise] + _b[out:exw]*(_b[med:_cons]
> + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5) + ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5))
        cde0: _b[out:exercise]*2.5
        pnie: _b[out:wtchg]*_b[med:exercise]*2.5
        tnde: (_b[out:exercise] + _b[out:exw]*((_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)) + _b[med:exe
> rcise]*2.5))*2.5

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         nde |  -.3191431    .792894    -0.40   0.687    -1.873187    1.234901
         nie |  -4.939743   .3844027   -12.85   0.000    -5.693158   -4.186327
          te |  -5.258886   .6638817    -7.92   0.000     -6.56007   -3.957701
          pm |   .9393135   .1441129     6.52   0.000     .6568574     1.22177
        cde0 |  -2.087915   .6905018    -3.02   0.002    -3.441273   -.7345559
        pnie |  -.8655719   .4624469    -1.87   0.061    -1.771951    .0408073
        tnde |  -4.393314   .6362454    -6.91   0.000    -5.640332   -3.146296
------------------------------------------------------------------------------

R: สูตรในรูปปิดชุดเดียวกันด้วยวิธีเดลตาเชิงตัวเลข

# Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
effects <- function(t) {
  m0 <- t[["med_(Intercept)"]] + t[["med_age"]] * agebar + t[["med_male"]] * malebar
  b1 <- t[["med_exercise"]]
  nde <- (t[["out_exercise"]] + t[["out_exw"]] * m0) * 2.5
  nie <- (t[["out_wtchg"]] + t[["out_exw"]] * 2.5) * b1 * 2.5
  c(nde = nde, nie = nie, te = nde + nie, pm = nie / (nde + nie),
    cde0 = t[["out_exercise"]] * 2.5,
    pnie = t[["out_wtchg"]] * b1 * 2.5,
    tnde = (t[["out_exercise"]] + t[["out_exw"]] * (m0 + b1 * 2.5)) * 2.5)
}
# Delta method with a central-difference Jacobian
est <- effects(th)
J <- sapply(seq_along(th), function(j) {
  h <- 1e-6 * max(1, abs(th[j]))
  up <- th; dn <- th
  up[j] <- up[j] + h
  dn[j] <- dn[j] - h
  (effects(up) - effects(dn)) / (2 * h)
})
se <- sqrt(diag(J %*% V %*% t(J)))
for (k in c("nde", "nie", "te", "pm", "cde0")) {
  canon(k, est[[k]])
  canon(paste0(k, ".lb"), est[[k]] - z975 * se[[k]])
  canon(paste0(k, ".ub"), est[[k]] + z975 * se[[k]])
}
canon("pnie", est[["pnie"]])
canon("tnde", est[["tnde"]])

> effects <- function(t) {
+     m0 <- t[["med_(Intercept)"]] + t[["med_age"]] * agebar +
+         t[["med_male"]] * malebar
+     b1 <- t[["med_exercise"]]
+     nde <- (t[["out_exercise"]] + t[["out_exw"]] * m0) * 2.5
+     nie <- (t[["out_wtchg"]] + t[["out_exw"]] * 2.5) * b1 * 2.5
+     c(nde = nde, nie = nie, te = nde + nie, pm = nie/(nde + nie),
+         cde0 = t[["out_exercise"]] * 2.5, pnie = t[["out_wtchg"]] *
+             b1 * 2.5, tnde = (t[["out_exercise"]] + t[["out_exw"]] *
+             (m0 + b1 * 2.5)) * 2.5)
+ }

> est <- effects(th)

> J <- sapply(seq_along(th), function(j) {
+     h <- 1e-06 * max(1, abs(th[j]))
+     up <- th
+     dn <- th
+     up[j] <- up[j] + h
+     dn[j] <- dn[j] - h
+     (effects(up) - effects(dn))/(2 * h)
+ })

> se <- sqrt(diag(J %*% V %*% t(J)))

> for (k in c("nde", "nie", "te", "pm", "cde0")) {
+     canon(k, est[[k]])
+     canon(paste0(k, ".lb"), est[[k]] - z975 * se[[k]])
+     canon(paste0(k, ".ub"), est[[k]] + z975 * se[[k]])
+ }
CANON w6.nde -0.3191
CANON w6.nde.lb -1.8732
CANON w6.nde.ub 1.2349
CANON w6.nie -4.9397
CANON w6.nie.lb -5.6932
CANON w6.nie.ub -4.1863
CANON w6.te -5.2589
CANON w6.te.lb -6.5601
CANON w6.te.ub -3.9577
CANON w6.pm 0.9393
CANON w6.pm.lb 0.6569
CANON w6.pm.ub 1.2218
CANON w6.cde0 -2.0879
CANON w6.cde0.lb -3.4413
CANON w6.cde0.ub -0.7346

> canon("pnie", est[["pnie"]])
CANON w6.pnie -0.8656

> canon("tnde", est[["tnde"]])
CANON w6.tnde -4.3933

Stata: คำสั่ง mediate ที่มีในตัว

capture noisily mediate (sbp age male) (wtchg age male) (exercise, continuous(0 2.5)), nie nde pnie tnde te

. capture noisily mediate (sbp age male) (wtchg age male) (exercise, continuous(0 2.5)), nie nde pnie tnde te

Iteration 0:  EE criterion =  2.075e-26
Iteration 1:  EE criterion =  9.065e-28

Causal mediation analysis                                Number of obs = 2,000

Outcome model:     Linear
Mediator model:    Linear
Mediator variable: wtchg
Treatment type:    Continuous

Continuous treatment levels:
  0: exercise =   0 (control)
  1: exercise = 2.5

------------------------------------------------------------------------------------
                   |               Robust
               sbp | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
NIE                |
          exercise |
         (1 vs 0)  |  -4.939743   .3481939   -14.19   0.000     -5.62219   -4.257295
-------------------+----------------------------------------------------------------
NDE                |
          exercise |
         (1 vs 0)  |  -.3191431   .7607765    -0.42   0.675    -1.810238    1.171951
-------------------+----------------------------------------------------------------
PNIE               |
          exercise |
         (1 vs 0)  |  -.8655719    .481407    -1.80   0.072    -1.809112    .0779685
-------------------+----------------------------------------------------------------
TNDE               |
          exercise |
         (1 vs 0)  |  -4.393314   .5850942    -7.51   0.000    -5.540077    -3.24655
-------------------+----------------------------------------------------------------
TE                 |
          exercise |
         (1 vs 0)  |  -5.258886   .6896579    -7.63   0.000     -6.61059   -3.907181
------------------------------------------------------------------------------------

Stata: คำสั่ง paramed ที่ผู้ใช้เขียนขึ้น

    capture noisily paramed sbp, avar(exercise) mvar(wtchg) cvars(age male) a0(0) a1(2.5) m(0) yreg(linear) mreg(linear)

.     capture noisily paramed sbp, avar(exercise) mvar(wtchg) cvars(age male) a0(0) a1(2.5) m(0) yreg(linear) mreg(linea
> r)

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(5, 1994)      =    358.74
       Model |  278452.465         5  55690.4931   Prob > F        =    0.0000
    Residual |  309547.372     1,994  155.239404   R-squared       =    0.4736
-------------+----------------------------------   Adj R-squared   =    0.4722
       Total |  587999.837     1,999  294.146992   Root MSE        =     12.46

-----------------------------------------------------------------------------------
              sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
------------------+----------------------------------------------------------------
         exercise |  -.8351658   .2762007    -3.02   0.003    -1.376838   -.2934936
            wtchg |   .2898677   .1545066     1.88   0.061    -.0131437    .5928791
_exercise_X_wtchg |   .5457527   .0547717     9.96   0.000     .4383369    .6531685
              age |   .4668235   .0221736    21.05   0.000     .4233377    .5103092
             male |   3.880244   .5609906     6.92   0.000     2.780054    4.980433
            _cons |   101.6947   1.502064    67.70   0.000     98.74896    104.6405
-----------------------------------------------------------------------------------

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(3, 1996)      =    412.53
       Model |  7937.34882         3  2645.78294   Prob > F        =    0.0000
    Residual |  12801.4031     1,996  6.41352861   R-squared       =    0.3827
-------------+----------------------------------   Adj R-squared   =    0.3818
       Total |  20738.7519     1,999  10.3745632   Root MSE        =    2.5325

------------------------------------------------------------------------------
       wtchg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.194437   .0435007   -27.46   0.000    -1.279749   -1.109126
         age |   .0236362   .0044603     5.30   0.000     .0148888    .0323836
        male |   .2156961   .1139233     1.89   0.058     -.007725    .4391173
       _cons |  -.1197085   .3035733    -0.39   0.693    -.7150622    .4756452
------------------------------------------------------------------------------

             |   Estimate     Std Err   P>|z|   [95% Conf   Interval]
-------------+-------------------------------------------------------
         cde | -2.0879146   .69050183   0.002  -3.4412982    -.734531
         nde | -.31914317   .79289403   0.687  -1.8732155   1.2349291
         nie | -4.9397425   .38440274   0.000  -5.6931719  -4.1863132
         mte | -5.2588857    .6638817   0.000  -6.5600938  -3.9576776

cde:controlled direct effect, nde:natural direct effect, nie:natural indirect effect, mte:marginal total effect

R: แพ็กเกจ mediation

# Cross-check with the mediation package (quasi-Bayesian simulation, 1000 draws)
set.seed(202664)
fy_int <- lm(sbp ~ exercise * wtchg + age + male, data = d)
med <- mediation::mediate(fm, fy_int, treat = "exercise", mediator = "wtchg",
                          control.value = 0, treat.value = 2.5, sims = 1000)
# d1 = indirect effect with exposure at a (NIE); z0 = direct effect with M at M(a*) (NDE)
canon("mediation.nie.r", med$d1)
canon("mediation.nie.lb.r", med$d1.ci[1])
canon("mediation.nie.ub.r", med$d1.ci[2])
canon("mediation.nde.r", med$z0)
canon("mediation.nde.lb.r", med$z0.ci[1])
canon("mediation.nde.ub.r", med$z0.ci[2])
canon("mediation.pnie.r", med$d0)
canon("mediation.tnde.r", med$z1)
canon("mediation.te.r", med$tau.coef)
canon("mediation.te.lb.r", med$tau.ci[1])
canon("mediation.te.ub.r", med$tau.ci[2])

> set.seed(202664)

> fy_int <- lm(sbp ~ exercise * wtchg + age + male,
+     data = d)

> med <- mediation::mediate(fm, fy_int, treat = "exercise",
+     mediator = "wtchg", control.value = 0, treat.value = 2.5,
+     sims = 1000)

> canon("mediation.nie.r", med$d1)
CANON w6.mediation.nie.r -4.9542

> canon("mediation.nie.lb.r", med$d1.ci[1])
CANON w6.mediation.nie.lb.r -5.7222

> canon("mediation.nie.ub.r", med$d1.ci[2])
CANON w6.mediation.nie.ub.r -4.1960

> canon("mediation.nde.r", med$z0)
CANON w6.mediation.nde.r -0.3032

> canon("mediation.nde.lb.r", med$z0.ci[1])
CANON w6.mediation.nde.lb.r -1.8452

> canon("mediation.nde.ub.r", med$z0.ci[2])
CANON w6.mediation.nde.ub.r 1.3059

> canon("mediation.pnie.r", med$d0)
CANON w6.mediation.pnie.r -0.8804

> canon("mediation.tnde.r", med$z1)
CANON w6.mediation.tnde.r -4.3771

> canon("mediation.te.r", med$tau.coef)
CANON w6.mediation.te.r -5.2574

> canon("mediation.te.lb.r", med$tau.ci[1])
CANON w6.mediation.te.lb.r -6.5889

> canon("mediation.te.ub.r", med$tau.ci[2])
CANON w6.mediation.te.ub.r -3.9108

ผลลัพธ์แบบใช่หรือไม่ใช่: คำถามเรื่องมาตรวัดกลับมาอีกครั้ง

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

บนมาตรวัดผลต่างความเสี่ยง (risk-difference scale) ผลตามธรรมชาติแต่ละตัวคือผลต่างระหว่างความเสี่ยงเฉลี่ยสองค่า เมื่อแบบจำลองของผลลัพธ์เป็นโลจิสติกจะไม่มีสูตรในรูปปิดที่เรียบง่าย ผลจึงคำนวณด้วยmediational g-formula สำหรับแต่ละคน แบบจำลองของผลลัพธ์ที่ประมาณแล้วทำนายความเสี่ยงของความดันโลหิตสูง โดยกำหนดการออกกำลังกายไว้ที่ $a^*$ หรือ $a$ และให้การเปลี่ยนแปลงของน้ำหนักกระจายตามการแจกแจงที่ประมาณได้ภายใต้ $a^*$ หรือ $a$ จากนั้นเฉลี่ยความเสี่ยงเหล่านี้ข้ามทุกคน ผลทั้งสามตัวเปรียบเทียบความเสี่ยงเฉลี่ยเป็นคู่ ๆ ตรงตามนิยามข้างต้น

บนมาตรวัด odds ratio (odds-ratio scale) มีสูตรในรูปปิดสำหรับผลลัพธ์แบบโลจิสติกที่ใช้ตัวแปรสื่อกลางเชิงเส้น [2] สูตรเหล่านั้นอาศัยการประมาณเมื่อผลลัพธ์พบน้อย (rare-outcome approximation) คือถือว่า odds ratio เป็นเหมือนอัตราส่วนความเสี่ยง ซึ่งใกล้เคียงก็ต่อเมื่อผลลัพธ์พบไม่บ่อยในทุกระดับของสิ่งสัมผัส ตัวแปรสื่อกลาง และตัวแปรร่วม ความดันโลหิตสูงใน cohort นี้ไม่ได้พบน้อย

ความดันโลหิตสูง: ค่าจริงในข้อมูลจำลอง ค่าประมาณ และช่วงเชื่อมั่นจากวิธีเดลตา

อ่านผลบนมาตรวัดทั้งสอง

บนมาตรวัดผลต่างความเสี่ยง ช่วงเชื่อมั่นจากวิธีเดลตาของผลทางตรงตามธรรมชาติ ผลทางอ้อมตามธรรมชาติ และผลรวมทั้งหมด ครอบคลุมค่าจริงในข้อมูลจำลองของแต่ละตัว ค่าประมาณของผลทางตรงตามธรรมชาติคือ -0.068 ซึ่งอยู่ไกลจาก -0.015 อีกครั้งเพราะส่วนที่เป็นผลทางตรงมีความแปรผันมากกว่า ช่วงเชื่อมั่นจากวิธีเดลตาของสัดส่วนของผลที่ผ่านตัวแปรสื่อกลาง คือ 0.31 ถึง 0.87 ไม่ครอบคลุมค่าจริงในข้อมูลจำลองที่ 0.89 ช่วงเชื่อมั่น 95% คาดว่าจะพลาดบ้างเป็นครั้งคราว และอัตราส่วนของผลต่างที่แปรผันสองตัวคือจุดที่การประมาณด้วยเส้นตรงของวิธีเดลตาอ่อนที่สุด

บนมาตรวัด odds ratio สูตรในรูปปิดเมื่อผลลัพธ์พบน้อยให้ 1.03 สำหรับผลทางตรงตามธรรมชาติ และ 0.51 สำหรับผลทางอ้อมตามธรรมชาติในประชากรจำลอง odds ratio ที่แท้จริงในประชากรนั้นที่สัดส่วนอายุและเพศเฉลี่ยคือ 0.93 และ 0.53 เมื่อผลลัพธ์พบบ่อยเช่นนี้ การประมาณยังทำให้ผลทางตรงไปอยู่คนละฝั่งของ 1

ค่าประมาณจากตัวอย่าง คือ 0.76 (ช่วงเชื่อมั่นจากวิธีเดลตา 0.53 ถึง 1.08) และ 0.58 (ช่วงเชื่อมั่นจากวิธีเดลตา 0.49 ถึง 0.68) ประมาณปริมาณที่เป็นค่าประมาณ ไม่ใช่ค่าที่แท้จริง แต่ละช่วงครอบคลุมเป้าหมายตามสูตรเมื่อผลลัพธ์พบน้อยของมัน คือ 1.03 และ 0.51 ที่ขนาดตัวอย่างนี้ช่วงเหล่านั้นยังครอบคลุมค่า exact ด้วย ข้อมูลเพียงอย่างเดียวจึงไม่แสดงว่าสูตรมุ่งไปที่ปริมาณที่ผิด

Stata: แบบจำลองโลจิสติกของผลลัพธ์ที่มีปฏิกิริยาสัมพันธ์

* Outcome model: logistic regression with the exposure-mediator interaction
logit htn exercise wtchg exw age male

. * Outcome model: logistic regression with the exposure-mediator interaction
. logit htn exercise wtchg exw age male

Iteration 0:  Log likelihood = -1241.3831
Iteration 1:  Log likelihood = -1015.5316
Iteration 2:  Log likelihood = -1000.2137
Iteration 3:  Log likelihood = -999.57991
Iteration 4:  Log likelihood = -999.57787
Iteration 5:  Log likelihood = -999.57787

Logistic regression                                     Number of obs =  2,000
                                                        LR chi2(5)    = 483.61
                                                        Prob > chi2   = 0.0000
Log likelihood = -999.57787                             Pseudo R2     = 0.1948

------------------------------------------------------------------------------
         htn | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -.1819008   .0538403    -3.38   0.001    -.2874259   -.0763757
       wtchg |   .1060281   .0314984     3.37   0.001     .0442925    .1677638
         exw |   .0314067   .0156946     2.00   0.045     .0006459    .0621675
         age |   .0514611   .0045331    11.35   0.000     .0425763    .0603459
        male |   .1396875   .1104594     1.26   0.206    -.0768089    .3561839
       _cons |  -3.431735   .3109377   -11.04   0.000    -4.041162   -2.822308
------------------------------------------------------------------------------

R: ผลต่างความเสี่ยงด้วย mediational g-formula

# Mediational g-formula on the risk-difference scale: average over people and 20 normal quantiles of M
zq <- qnorm((seq_len(20) - 0.5) / 20)
i <- rep(seq_len(nrow(db)), each = 20)
z <- rep(zq, times = nrow(db))
m0 <- cm[["(Intercept)"]] + cm[["age"]] * db$age[i] + cm[["male"]] * db$male[i] + sig * z
m1 <- m0 + cm[["exercise"]] * 2.5
xc <- g[["(Intercept)"]] + g[["age"]] * db$age[i] + g[["male"]] * db$male[i]
r00 <- mean(plogis(xc + g[["wtchg"]] * m0))
r10 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] * m0 + g[["exw"]] * 2.5 * m0))
r11 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] * m1 + g[["exw"]] * 2.5 * m1))
canon("bin.rd_nde", r10 - r00)
canon("bin.rd_nie", r11 - r10)
canon("bin.rd_te", r11 - r00)
canon("bin.rd_pm", (r11 - r10) / (r11 - r00))
# delta-method 95% CIs on the risk-difference scale (NDE, NIE, TE, proportion mediated)
for (k in 1:4) {
  e <- c("rd_nde", "rd_nie", "rd_te", "rd_pm")[k]
  canon(paste0("bin.", e, ".lb"), f0b[k] - z975 * seb[k])
  canon(paste0("bin.", e, ".ub"), f0b[k] + z975 * seb[k])
}

> zq <- qnorm((seq_len(20) - 0.5)/20)

> i <- rep(seq_len(nrow(db)), each = 20)

> z <- rep(zq, times = nrow(db))

> m0 <- cm[["(Intercept)"]] + cm[["age"]] * db$age[i] +
+     cm[["male"]] * db$male[i] + sig * z

> m1 <- m0 + cm[["exercise"]] * 2.5

> xc <- g[["(Intercept)"]] + g[["age"]] * db$age[i] +
+     g[["male"]] * db$male[i]

> r00 <- mean(plogis(xc + g[["wtchg"]] * m0))

> r10 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] *
+     m0 + g[["exw"]] * 2.5 * m0))

> r11 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] *
+     m1 + g[["exw"]] * 2.5 * m1))

> canon("bin.rd_nde", r10 - r00)
CANON w6.bin.rd_nde -0.0679

> canon("bin.rd_nie", r11 - r10)
CANON w6.bin.rd_nie -0.0988

> canon("bin.rd_te", r11 - r00)
CANON w6.bin.rd_te -0.1667

> canon("bin.rd_pm", (r11 - r10)/(r11 - r00))
CANON w6.bin.rd_pm 0.5925

> for (k in 1:4) {
+     e <- c("rd_nde", "rd_nie", "rd_te", "rd_pm")[k]
+     canon(paste0("bin.", e, ".lb"), f0b[k] - z975 * seb[k])
+     canon(paste0("bin.", e, ".ub"), f0b[k] + z975 * seb[k])
+ }
CANON w6.bin.rd_nde.lb -0.1317
CANON w6.bin.rd_nde.ub -0.0042
CANON w6.bin.rd_nie.lb -0.1300
CANON w6.bin.rd_nie.ub -0.0676
CANON w6.bin.rd_te.lb -0.2170
CANON w6.bin.rd_te.ub -0.1164
CANON w6.bin.rd_pm.lb 0.3116
CANON w6.bin.rd_pm.ub 0.8735

สิ่งที่ต้องมีจึงจะอ่านผลเหล่านี้เชิงเหตุผลได้

นิยามข้างต้นเป็นเชิงเหตุผล ส่วนการถดถอยไม่ใช่ การก้าวจากอย่างหนึ่งไปอีกอย่างหนึ่งต้องมีข้อสมมติสี่ข้อ โดย $C$ คือตัวแปรร่วมที่วัดได้ [2, 5]

1. ไม่มีตัวกวนที่ไม่ได้วัดระหว่างสิ่งสัมผัสกับผลลัพธ์ เมื่อกำหนด $C$
2. ไม่มีตัวกวนที่ไม่ได้วัดระหว่างตัวแปรสื่อกลางกับผลลัพธ์ เมื่อกำหนด $A$ และ $C$
3. ไม่มีตัวกวนที่ไม่ได้วัดระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง เมื่อกำหนด $C$
4. ข้อสมมติข้ามโลก (cross-world assumption) เมื่อกำหนด $C$ ผลลัพธ์ที่อาจเกิดขึ้น $Y^{a,m}$ เป็นอิสระจากตัวแปรสื่อกลางที่อาจเกิดขึ้น $M^{a^*}$ ภายใต้อีกระดับของสิ่งสัมผัส เขียนเป็น $Y^{a,m} \perp M^{a^*} \mid C$

ผลทางตรงแบบควบคุมต้องการเพียงสองข้อแรก ผลตามธรรมชาติต้องการทั้งสี่ข้อ เพราะมี $Y^{a,M^{a^*}}$ ซึ่งเชื่อมการออกกำลังกายที่ $a$ เข้ากับการเปลี่ยนแปลงของน้ำหนักจากโลกที่การออกกำลังกายเป็น $a^*$

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

ทำไมจึงตรวจข้อสมมติข้ามโลกไม่ได้

แม้ไม่มีตัวแปรเช่นนั้น ก็ยังตรวจข้อสมมติข้ามโลกจากข้อมูลไม่ได้ ไม่มีใครถูกสังเกตได้ว่าออกกำลังกายที่ $a$ และในเวลาเดียวกันก็มีการเปลี่ยนแปลงของน้ำหนักเท่ากับที่เขาจะมีที่ $a^*$ การทดลองใด ๆ ก็สร้างการจับคู่เช่นนั้นไม่ได้เช่นกัน การสุ่มการออกกำลังกายรับประกันข้อสมมติข้อแรกและข้อที่สาม แต่ไม่รับประกันข้อที่สองและข้อที่สี่

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

การวิเคราะห์ความไวสำหรับตัวกวนที่ไม่ได้วัดระหว่างตัวแปรสื่อกลางกับผลลัพธ์

การสุ่มสิ่งสัมผัสรับประกันข้อสมมติข้อแรกและข้อที่สาม แต่ปล่อยข้อที่สองและข้อที่สี่ไว้เปิด ผู้ที่น้ำหนักลดมากกว่าอาจต่างจากผู้อื่นในด้านอื่นที่ลดความดันได้เช่นกัน เช่นการนอนหรือการดื่มแอลกอฮอล์ การวิเคราะห์ความไว (sensitivity analysis) ถามว่าตัวกวนที่ไม่ได้วัดเช่นนั้น คือ $U$ ต้องแรงเพียงใดจึงจะเปลี่ยนข้อสรุปได้

มีสองวิธีที่ใช้กันทั่วไป สูตรอคติ (bias formulas) [5] ให้นักวิเคราะห์ระบุว่า $U$ ส่งผลต่อความดันโลหิตแรงเพียงใด และการแจกแจงของ $U$ ต่างกันอย่างไรระหว่างระดับการออกกำลังกายในผู้ที่มีการเปลี่ยนแปลงของน้ำหนักเท่ากัน แล้วสูตรจะคืนผลทางตรงและทางอ้อมตามธรรมชาติที่แก้สำหรับตัวกวนระดับนั้น วิธีของ Imai และคณะ [6] ใช้สหสัมพันธ์ระหว่างพจน์คลาดเคลื่อนของแบบจำลองของตัวแปรสื่อกลางและของผลลัพธ์ ซึ่งเป็นศูนย์เมื่อไม่มีตัวกวนเช่นนั้น และรายงานว่าสหสัมพันธ์นั้นต้องมากเพียงใดผลทางอ้อมจึงจะหายไป

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

ขั้นตอนจากคำถามถึงรายงาน

ขั้นตอนด้านล่างเรียงนิยาม แบบจำลอง และข้อสมมติตามลำดับที่โปรโตคอลต้องการ

1. วาด DAG วางสิ่งสัมผัส ตัวแปรสื่อกลาง และผลลัพธ์ เพิ่มตัวกวนของความสัมพันธ์ทั้งสามคู่ และทำเครื่องหมายตัวกวนระหว่างตัวแปรสื่อกลางกับผลลัพธ์ที่สิ่งสัมผัสส่งผลกระทบ
2. ระบุ estimand คือปริมาณที่แน่ชัดที่ต้องการประมาณ ได้แก่ผล (ควบคุมที่ $m$ ที่ระบุ หรือตามธรรมชาติ) ระดับสิ่งสัมผัสสองระดับ $a^*$ และ $a$ มาตรวัด และประชากร
3. ประมาณแบบจำลองที่มีปฏิกิริยาสัมพันธ์ คงพจน์ผลคูณระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลางไว้ในแบบจำลองของผลลัพธ์ เว้นแต่ได้ระบุเหตุผลที่จะตัดออกไว้ล่วงหน้า
4. คำนวณผล ใช้สูตรในรูปปิดเมื่อแบบจำลองทั้งสองเป็นเชิงเส้น และใช้สูตร odds ratio เฉพาะเมื่อผลลัพธ์พบน้อย นอกนั้นใช้ mediational g-formula หรือการจำลอง
5. ใส่ช่วงเชื่อมั่น ช่วงเชื่อมั่นจากวิธีเดลตาทำได้เร็ว ช่วงเชื่อมั่นจากบูตสแตรป (bootstrap) ซึ่งประมาณแบบจำลองใหม่จากสำเนาข้อมูลที่สุ่มซ้ำจำนวนมาก เป็นอีกทางเลือกที่ใช้บ่อย โดยเฉพาะกับอัตราส่วน เช่นสัดส่วนของผลที่ผ่านตัวแปรสื่อกลาง
6. ทำการวิเคราะห์ความไว สำหรับตัวกวนที่ไม่ได้วัดระหว่างตัวแปรสื่อกลางกับผลลัพธ์
7. รายงาน DAG, estimand, แบบจำลองทั้งสอง, ผลพร้อมช่วงเชื่อมั่น และข้อสมมติ ตามแนวทางการรายงาน AGReMA สำหรับการวิเคราะห์ตัวแปรสื่อกลาง [8]

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

สัดส่วนของผลที่ผ่านตัวแปรสื่อกลาง ต้องอ่านด้วยความระมัดระวัง

สัดส่วนของผลที่ผ่านตัวแปรสื่อกลางคือ $\mathrm{NIE}/\mathrm{TE}$ บนมาตรวัดเดียว ในตัวอย่างคำนวณด้วยมือ ค่านี้คือ 0.83 จากการลดลงทั้งหมดของ SBP 5.95 mmHg มี 4.95 mmHg ที่ผ่านน้ำหนักในการแยกแบบผลตามธรรมชาติ มีคุณสมบัติสามข้อที่จำกัดว่าตัวเลขนี้อ่านได้ไกลเพียงใด

ประการแรก เมื่อมีปฏิกิริยาสัมพันธ์ระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลาง การแยกผลไม่ได้มีทางเดียว [3] ผลรวมทั้งหมดยังเท่ากับผลทางอ้อมแบบ pure (pure indirect effect, PNIE) ซึ่งคงการออกกำลังกายไว้ที่ $a^*$ ขณะที่น้ำหนักเลื่อน บวกผลทางตรงแบบ total (total direct effect, TNDE) ซึ่งคงการเปลี่ยนแปลงของน้ำหนักไว้ที่ระดับภายใต้ $a$ ตามภาษานี้ ผลทางตรงและทางอ้อมตามธรรมชาติที่ใช้มาจนถึงตอนนี้คือผลทางตรงแบบ pure และผลทางอ้อมแบบ total

ในที่นี้ เมื่อ $a^* = 0$ จะได้ $\mathrm{PNIE} = \theta_2 \beta_1 (a - a^*)$ เท่ากับ $0.4 \times (-3.0) = -1.20$ mmHg ซึ่งใช้ความชันของน้ำหนักเมื่อไม่ออกกำลังกาย ผลทางตรงแบบ total ใช้ค่าเฉลี่ยการเปลี่ยนแปลงของน้ำหนักภายใต้ $a$ ซึ่งคือ $1.2 + (-1.2)(2.5) = -1.8$ กิโลกรัม ความชันของการออกกำลังกายของมันจึงเป็น $-1 + 0.5 \times (-1.8) = -1.9$ และ $\mathrm{TNDE} = -1.9 \times 2.5 = -4.75$ mmHg เมื่อวัดด้วยผลทางอ้อมแบบ pure สัดส่วนของผลที่ผ่านตัวแปรสื่อกลางคือ $-1.20/-5.95 = 0.20$ ไม่ใช่ 0.83 จากแบบจำลองเดียวกันและคนกลุ่มเดียวกัน

ประการที่สอง สัดส่วนของผลที่ผ่านตัวแปรสื่อกลางเป็นอัตราส่วน และอัตราส่วนไม่เสถียรเมื่อตัวหารเล็กหรือไม่แน่นอน ใน cohort จำลอง ค่าประมาณคือ 0.94 โดยช่วงเชื่อมั่นจากวิธีเดลตาคือ 0.66 ถึง 1.22 ซึ่งครอบคลุมค่าจริงในข้อมูลจำลองที่ 0.83 และเลยไปเกิน 1 สัดส่วนที่เกิน 1 จะหมายความว่าผลทางตรงชี้ไปอีกทิศ ช่วงนี้จึงครอบคลุมเรื่องราวที่ต่างกันมาก

ประการที่สาม สัดส่วนของผลที่ผ่านตัวแปรสื่อกลางไม่ได้บอกว่าจะเกิดอะไรขึ้นหากตัดตัวแปรสื่อกลางออก คำถามนั้นตรึงการเปลี่ยนแปลงของน้ำหนักด้วยการแทรกแซง และผลทางตรงแบบควบคุมคือสิ่งที่ตอบ เมื่อตรึงการเปลี่ยนแปลงของน้ำหนักไว้ที่ 0 กิโลกรัมสำหรับทุกคน การออกกำลังกายที่เพิ่มขึ้นจะลด SBP ลง 2.50 mmHg ในประชากรจำลอง เทียบกับผลรวมทั้งหมด 5.95 mmHg

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

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

• เขียนคอนทราสต์ ($a^*$ และ $a$) ผล (แบบควบคุมหรือตามธรรมชาติ) และมาตรวัดไว้ในแผนการวิเคราะห์ก่อนประมาณแบบจำลองใด ๆ
• คงพจน์ผลคูณระหว่างสิ่งสัมผัสกับตัวแปรสื่อกลางไว้ในแบบจำลองของผลลัพธ์ เว้นแต่มีเหตุผลที่ระบุไว้ล่วงหน้าให้ตัดออก และรายงานค่าประมาณของพจน์นั้น
• รายงานผลทางตรงตามธรรมชาติ ผลทางอ้อมตามธรรมชาติ และผลรวมทั้งหมดไปด้วยกัน แต่ละผลพร้อมช่วงเชื่อมั่นและวิธีที่อยู่เบื้องหลังช่วงเชื่อมั่นนั้น
• สำหรับผลลัพธ์แบบสองค่าที่พบบ่อย ให้พิจารณามาตรวัดผลต่างความเสี่ยงร่วมกับ mediational g-formula แทนสูตร odds ratio เมื่อผลลัพธ์พบน้อย
• ระบุรายการตัวกวนระหว่างตัวแปรสื่อกลางกับผลลัพธ์ที่ปรับแล้ว ระบุตัวที่สิ่งสัมผัสส่งผลต่อ และวางแผนการวิเคราะห์ความไวสำหรับตัวที่วัดไม่ได้
• พิจารณารายงานผลทางตรงแบบควบคุมด้วย เมื่อตัวแปรสื่อกลางเป็นสิ่งที่กำหนดได้ด้วยการแทรกแซง เช่นขนาดยา