两因素方差分析 (Two-way ANOVA) 一(无交互作用)——Stata软件实现

关键词：Stata; 两因素方差分析; 交互作用; 主效应; 单独效应; 简单效应

一、案例介绍

研究A、B两种镇痛药物联合运用在癌症患者的镇痛效果。药物A取3个剂量：1.0、2.5、5.0，对应变量drug_a的1、2、3组；药物B也取3个剂量：5、15、30 μg，对应变量drug_b的1、2、3组，共9个处理组。将27名研究对象随机分成9组，每组3名，记录每名对象的镇痛时间Time (min)。试分析A、B两药联合运用的镇痛效果。部分数据见图1。本文案例可从“附件下载”处下载。

二、问题分析

本案例的分析目的是分析A、B两药联合运用的镇痛效果。临床上，药物之间联合运用往往会相互影响，这种影响可能为正向的增强效应，也可能为反向的拮抗作用。针对这种情况，可以使用多因素方差分析。由于本案例为药物A和药物B两个因素，因此可以使用两因素方差分析。但需要满足6个条件：

条件1：观察变量唯一，且为连续变量。本研究中观察变量只有镇痛时间，且为连续变量，该条件满足。

条件2：有两个因素，且都为分类变量。本研究中有药物A、药物B两个因素，都为分类变量，该条件满足。

条件3：观测值相互独立。本研究中各研究对象的观测值都是独立的，不存在互相干扰的情况，该条件满足。

条件4：观察变量不存在显著的异常值，该条件需要通过软件分析后判断。

条件5：各组、各水平观测值为正态(或近似正态)分布，该条件需要通过软件分析后判断。

条件6：相互比较的各处理水平(组别)的总体方差相等，即方差齐同，该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 适用条件判断

两因素方差分析过程中使用残差判断异常值，先使用两因素方差分析生成因变量残差。

1. 生成因变量残差

(1) 软件操作

anova time drug_a drug_b drug_a#drug_b
predict e, residual 

(2) 结果解读

通过“Data Editor”视图可以看到生成了一列新的变量“e (残差)”(图2)。

2. 条件4判断(异常值判断)

(1) 软件操作

绘制残差箱线图(图3)。

graph box e, over(drug_a, relabel(1 "1.0mg" 2 "2.5mg" 3 "5.0mg")) over(drug_b, relabel(1 "5μg" 2 "15μg" 3 "30μg"))

(2) 结果解读

图3残差的箱线图未提示任何异常值和极端值，满足条件4。

3. 条件5判断(正态性检验)

(1) 软件操作

两因素方差分析时，需要分别考察每一组数据的正态性或使用残差考察整体数据的正态性。

①分别考察每一组数据的正态性，结果如图4-1、图4-2所示。

bysort drug_a: swilk time
bysort drug_b: swilk time

②使用残差考察整体数据的正态性，结果如图5所示。

swilk e

(2) 结果解读

图4-1、图4-2的正态性检验结果显示各组的P值均>0.1，提示各组数据均服从正态分布。图5的残差正态性检验结果显示，P值>0.1，提示因变量残差整体服从正态分布。综上，本案例满足条件5。此外，也可以使用Q-Q图判断正态性(读者可自行操作)。

4. 条件6判断(方差齐性检验)

(1) 软件操作

①绘制残差图检验方差齐性，结果如图6所示。

rvfplot, yline(0)

②使用BP(Breusch-Pagan)法检验方差齐性，结果如图7所示。

estat hettest e drug_a drug_b

(2) 结果解读

由图6可见，残差均匀地分布在其均值的上下两侧，提示观察变量的残差满足方差齐性；根据图7结果，P>0.1，不能拒绝方差相等的原假设，提示观察变量的残差满足方差齐性。综上，本案例满足条件6。

(二) 统计描述

1. 软件操作

tabulate drug_a drug_b, summarize(time)

2. 结果解读

图8的结果中列出了各组的均值和标准差，可知药物A在1.0 mg水平(drug_a=1)，药物B为5、15、30 μg水平(drug_b=1、2、3)时，镇痛时间分别为：79.000±4.000、94.333±2.517、108.667±4.509 min；药物A在2.5 mg水平(drug_a=2)，药物B为5、15、30 μg水平(drug_b=1、2、3)时，镇痛时间分别为：90.667±4.163、106.667±4.041、116.667±2.517 min；药物A在5.0 mg水平(drug_a=3)，药物B为5、15、30 μg水平(drug_b=1、2、3)时，镇痛时间分别为：108.000±4.583、118.333±5.508、132.000±4.000 min。

(三) 两因素方差分析

1. 软件操作

margins drug_a, over(drug_b) plot

2. 结果解读

图9的结果为药物A和药物B估算边际均值的相关统计量结果；图10绘制了各组镇痛时间的变化情况，可见随着药物A浓度的增加，镇痛时间逐渐延长；随着药物B浓度的增加，镇痛时间也逐渐延长；并且增加的幅度基本保持一致。

(四) 交互作用判断

由于本案例有两个因素(一个为药物A，另一个为药物B)，因此需要首先判断两个因素(药物)之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。

1. 软件操作

anova time drug_a drug_b drug_a#drug_b

2. 结果解读

图11的结果显示，drug_a与drug_b之间的交互作用结果为F_{drug_a#drug_b}=0.70，P=0.603，提示药物A与药物B之间的交互作用无统计学意义；药物A的F_{drug_a}=88.16，P <0.001，认为不同药物浓度时，镇痛时间差异有统计学意义；药物B的F_{drug_b}=95.27，P <0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。因此，本案例可直接以主效应分析结果进行判断，如果交互作用有统计学意义，则需要分析单独效应。

(五) 事后检验(两两比较)

上面分析得出了药物A、药物B在不同药物浓度时，镇痛时间差异均有统计学意义的结论，因此需要进行事后检验，开展两两比较。

1. 软件操作

pwcompare drug_a, mcompare(bonferroni) effects
pwcompare drug_b, mcompare(bonferroni) effects

2. 结果解读

图12-1显示了药物A不同浓度间两两比较的相关统计量结果。可知，2.5 mg(drug_a=2)和5.0 mg(drug_a=3)与1.0 mg(drug_a=1)相比，均值差逐渐增大。2.5 mg比1.0 mg时均值增加了10.667 min；至5.0 mg时均值增加了25.444 min；2.5 mg到5.0 mg时均值增加了14.778 min；以上差异均有统计学意义(P<0.001)。表明：随着药物A药物浓度的增加，镇痛时间呈上升趋势。

图12-2显示了药物B不同浓度间两两比较的相关统计量结果。可知，15.0 μg(drug_b=2)和30.0 μg(drug_b=3)与5.0 μg(drug_b=1)相比，均值差逐渐增大。15.0 μg比5.0 μg时均值增加了13.889 min；至30.0 μg时均值增加了26.556 min；15.0 μg到30.0 μg时均值增加了12.667 min；以上差异均有统计学意义(P<0.001)。表明：随着药物B药物浓度的增加，镇痛时间呈上升趋势。

四、结论

本研究采用两因素方差分析探讨A、B两种镇痛药物联合运用在癌症患者的镇痛效果。通过对模型残差绘制箱线图提示，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过对残差绘制残差图和进行BP(Breusch-Pagan)检验，提示数据总体方差相等；方差分析结果显示两药物之间无交互作用(F_{drug_a#drug_b} =0.70，P=0.603)，故进行主效应分析。

药物A在1.0 mg水平，药物B为5、15、30 μg水平时，镇痛时间分别为：79.000±4.000、94.333±2.517、108.667±4.509 min；药物A在2.5 mg水平，药物B为5、15、30 μg水平时，镇痛时间分别为：90.667±4.163、106.667±4.041、116.667±2.517 min；药物A在5.0 mg水平，药物B为5、15、30 μg水平时，镇痛时间分别为：108.000±4.583、118.333±5.508、132.000±4.000 min。药物A的主效应检验结果为，F_{drug_a}=88.16，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。药物B的主效应检验结果为，F_{drug_b} =95.27，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。进一步采用“Bonferroni”校正法进行两两比较可知，药物A 2.5 mg和5.0 mg与1.0 mg相比，均值差逐渐增大，差异均有统计学意义(P<0.001)；表明：随着药物A药物浓度的增加，镇痛时间呈上升趋势。药物B15.0 μg和30.0 μg与5.0 μg相比，均值差逐渐增大，差异均有统计学意义(P<0.001)；表明：随着药物B药物浓度的增加，镇痛时间呈上升趋势。综上可知，不同剂量药物A的镇痛效果不同，不同剂量药物B的镇痛效果不同，A、B两种药物不存在交互作用。

五、分析小技巧

• 正态性检验：两因素或多因素方差分析时，有两种选择来测试正态性：如果每组有较多观察数，且组别较少时，可使用原始数据检查每个组的正态性。如果有很多组，或每个组的观察数很少，可使用残差检查整体的正态性。在Stata中使用“swilk”命令可进行“Shapiro-Wilk (夏皮罗-威尔克正态性)”检验。关于正态性检验的注意事项详见文章(医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现)。
• 方差齐性检验：两因素方差分析时，使用残差检验方差齐性。可绘制残差图(“rvfplot”命令)或使用BP(Breusch-Pagan)法(“estat hettest”命令)检验方差齐性。关于方差齐性检验的更多内容请阅读(医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现)。
• 交互作用判断：两因素方差分析时，需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。当不存在交互作用时，说明两因素的作用彼此独立，逐一分析各因素的主效应即可；计算主效应时，在模型中仍需要保留交互项。
• 在上面的分析中，drug_a和drug_b取值都是1、2、3，得到的结果也是以1、2、3显示。根据案例介绍可知，药物A取3个剂量：1.0 、2.5、5.0 mg，对应变量drug_a的1、2、3组；药物B也取3个剂量：5、15、30 μg，对应变量drug_b的1、2、3组。如果想要在Stata分析结果中直接显示剂量值，可使用“recode”命令。

recode drug_a (1=1 "1.0mg") (2=2 "2.5mg") (3=3 "5.0mg"), pre(new)
recode drug_b (1=1 "5μg") (2=2 "15μg") (3=3 "30μg"), pre(new)

上述命令生成两个新变量newdrug_a(取值为1.0、2.5、5.0mg)和newdrug_b(取值为5、15、30μg)，分别对应drug_a和drug_b，部分数据如图13所示。后续分析时使用newdrug_a和newdrug_b即可，读者可自行尝试。