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

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

一、案例介绍

研究A、B两种镇痛药物联合运用在癌症患者的镇痛效果。A药取3个剂量：1.0 mg、2.5 mg、5.0 mg；B药也取3个剂量：5 μg、15 μg、30 μg，共9个处理组。将27名研究对象随机分成9组，每组3名，记录每名对象的镇痛时间Time (min)。试分析A、B两药联合运用的镇痛效果。

创建代表处理因素的变量“Drug_A”和“Drug_B”， “Drug_A”赋值为“1”、“2”、“3”分别代表1.0 mg、2.5 mg、5.0 mg三个剂量，“Drug_B”赋值为“1”、“2”、“3”分别代表5 μg、15 μg、30 μg三个剂量。创建观察变量“Time”，记录各处理组中每名对象的镇痛时间。部分数据见图1。本文案例可从“附件下载”处下载。

二、问题分析

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

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

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

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

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

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

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

三、软件操作及结果解读

(一) 导入数据

①利用LIBNAME语句建立SAS逻辑库关联，注意逻辑库名称要求，即最大长度8字符，必须以字母或下划线“_”开始，可以是字母、数字和下划线的任意组合。具体代码如下：

libname mydata 'D:\mydata';

通过这一步骤，SAS能够识别引号中的物理位置，将逻辑库建立在该目录下，同时在以下过程中新建的SAS表格便可以永久储存在该位置，便于反复读取和使用。先运行该代码使其生效。

②利用PROC IMPORT语句导入文件，代码如下：

proc import out= mydata.example
datafile=" D:\mydata\两因素方差分析一.csv "
dbms=csv replace;
getnames=yes;
run;

该过程在mydata逻辑库中生成example数据集，数据文件由DATAFILE=选项指定，DBMS=选项指定其数据库类型。该案例中初始数据集为csv文件，故而使用“dbms=csv”指定。如果已经存在相同名称的SAS数据集，即可使用REPLACE选项进行覆盖。GERNAMES=YES选项指定从第2行开始读取数据，将数据集的首行变量名作为SAS数据集的变量名。

(二) 适用条件判断

1. 生成因变量残差

ods output MEANS=resid;
proc anova data = mydata.example;
	class Drug_A Drug_B;
   model time = Drug_A|Drug_B;
	means Drug_A|Drug_B;
run;
quit;

proc sql;
create table resid_1 as
select Drug_A,Drug_B,Mean_Time
from resid
where n=3
order by Drug_A,Drug_B;
quit;

proc sort data=mydata.example;
by Drug_A Drug_B;
proc print;run;

data resid_1;
set resid_1(rename=(Drug_A=Drug_A1 Drug_B=Drug_B1));
Drug_A=input(Drug_A1,8.);
Drug_B=input(Drug_B1,8.);
drop Drug_A1 Drug_B1;
proc print;run;

data mydata.example_1;
merge mydata.example resid_1;
by Drug_A Drug_B;
resid=Time-Mean_Time;
proc print;run;

生成新数据集mydata.example_1，包含残差变量resid。

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

(1) 软件操作

异常值借助箱线图判断，使用proc sgplot可以绘制每一组镇痛时间Time残差的箱线图，具体代码如下所示：

proc sgplot data=mydata.example_1;
vbox resid/category=Drug_B group=Drug_A;
run;

(2) 结果解读

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

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

(1) 软件操作

proc univariate data=mydata.example_1 normal;
	class Drug_A Drug_B;
	var resid;
run;

其中，normal选项表示对变量做正态性分析。

(2) 结果解读

图3显示了“Tests for Normality (正态性检验)”的结果，包括Shapiro-Wilk (夏皮罗-威尔克正态性，S-W)检验和Kolmogorov-Smirnov (柯尔莫哥洛夫-斯米诺夫，K-S)检验等。K-S检验适用于大样本资料，本案查看S-W检验结果，可见Drug_A 1.0 mg和Drug_B 5 μg的处理组的P值>0.1，提示所有该处理组数据服从正态分布。其他8个处理组的情况也类似(因篇幅原因，不做额外展示)，均服从正态分布。其他处理组的Q-Q图也显示散点基本围绕对角线分布，提示数据服从正态分布。综上，本案例满足条件5。

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

(1) 软件操作

proc anova data = mydata.example_1;
class Drug_A Drug_B;
model resid = Drug_A*Drug_B;
means Drug_A*Drug_B/hovtest=levene(type=abs);
run;
quit;

(2) 结果解读

图4 “Levene's Test for Homogeneity (Levene’s方差齐性检验)”的结果显示，统计量F=0.30、P=0.9546；提示数据总体方差相等。本案例满足条件6。

(二) 统计描述及推断

1. 软件操作

使用如下代码对每一组的镇痛时间Time进行统计描述：

proc means data = mydata.example;
class Drug_A Drug_B;
var time;
run;

使用如下代码进行统计推断、交互作用判断以及主效应分析：

proc anova data = mydata.example;
class Drug_A Drug_B;
model time = Drug_A|Drug_B;
run;

model语句中的“Drug_A|Drug_B“等价于“Drug_A Drug_B Drug_A*Drug_B。

2. 结果解读

(1) 统计描述

图5 “Analysis Variable : Time (分析变量：Time)”表格中列出了各组的均值和标准差，可知A药在1.0 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：79.00±4.00 min、94.33±2.52 min、108.67±4.51 min；A药在2.5 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：90.67±4.16 min、106.67±4.04 min、116.67±2.52 min；A药在5.0 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：108.00±4.58 min、118.33±5.51 min、132.00±4.00 min。

(2) 统计推断

图6显示了对所假设的模型进行方差分析的结果，原假设为模型中所有的影响因素均无作用，即A药、B药、两者的交互作用均对镇痛时间无影响。第一行的“Model(模型)”即对所假设的模型进行检验的结果，F=46.21，P<0.001，因此所用的模型有统计学意义，以上所提到的影响因素中至少有一个对镇痛时间有影响，但具体是哪些需看后续分析结果。

(3) 交互作用判断

由于本案例有两个因素(一个为A药，另一个为B药)，因此首先需要判断两个药物之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。

图7第三行是对A药和B药的交互作用的判断，P值为0.6027，无统计学意义，说明两者不存在交互作用。

(4) 主效应分析

图7第一行是对A药的主效应检验，结果显示F_{Drug_A}=88.16，P<0.001，认为A药不同浓度时，镇痛时间差异有统计学意义。

图7第二行是对B药的主效应检验，结果显示F_{Drug_B}=95.27，P<0.001，认为B药不同浓度时，镇痛时间差异有统计学意义。

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

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

1. 软件操作

运用GLM过程进行事后检验，具体代码如下所示：

proc glm data = mydata.example;
class Drug_A Drug_B;
model time = Drug_A Drug_B;
lsmeans Drug_A Drug_B/ tdiff adjust=bon;
run;
quit;

class语句声明分类变量。Model语句指明模型因变量和自变量，由于A药和B药不存在交互作用，因此，模型只保留Drug_A和Drug_B。 LSMEANS语句分别对A药和B药不同浓度间进行两两比较，tdiff选项输出两两比较的T统计值和P值，adjust选项指明使用Bonferroni方法进行两两比较。quit语句结束GLM过程。

2. 结果解读

图8提供了A药不同浓度镇痛时间Time的均值，以及不同浓度间两两比较的T统计值和P值。可知，2.5 mg和5.0 mg与1.0 mg相比，镇痛时间显著增长(T=5.71, P<0.0001; T=13.60, P<0.0001)。5.0 mg与2.5 mg相比，镇痛时间显著增长(T=7.90, P<0.0001)。该结果表明，随着A药药物浓度的增加，镇痛时间呈上升趋势。

图9提供了B药不同浓度镇痛时间Time的均值，以及不同浓度间两两比较的T统计值和P值。可知，15 μg和30 μg与5 μg相比，镇痛时间显著增长(T=7.42, P<0.0001; T=14.19, P<0.0001)。30 μg与15 μg相比，镇痛时间显著增长(T=6.77, P<0.0001)。该结果表明，随着B药药物浓度的增加，镇痛时间呈上升趋势。

四、结论

本研究采用两因素方差分析探讨A、B两种镇痛药物联合运用在癌症患者的镇痛效果。通过专业知识判断，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过Levene’s检验，提示数据总体方差相等；两药物之间无交互作用(F_{Drug_A*Drug_B}=0.70，P=0.6027)，故进行主效应分析。

A药在1.0 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：79.00±4.00 min、94.33±2.52 min、108.67±4.51 min；A药在2.5 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：90.67±4.16 min、106.67±4.04 min、116.67±2.52 min；A药在5.0 mg水平，B药为5 μg、15 μg、30 μg水平时，镇痛时间分别为：108.00±4.58 min、118.33±5.51 min、132.00±4.00 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相比，镇痛时间显著增长(T=5.71, P<0.0001; T=13.60, P<0.0001)；5.0 mg与2.5 mg相比，镇痛时间显著增长(T=7.90, P<0.0001)；表明随着A药药物浓度的增加，镇痛时间呈上升趋势。B药15 μg和30 μg与5 μg相比，镇痛时间显著增长(T=7.42, P<0.0001; T=14.19, P<0.0001)；30 μg与15 μg相比，镇痛时间显著增长(T=6.77, P<0.0001)；表明随着B药药物浓度的增加，镇痛时间呈上升趋势。综上可知，不同剂量药物A的镇痛效果不同，不同剂量药物B的镇痛效果不同，A、B两种药物不存在交互作用。

五、分析小技巧

• 正态性检验：两因素或多因素方差分析时，有两种选择来测试正态性：如果每组有较多观察数，且组别较少时，可使用原始数据检查每个组的正态性。如果有很多组，或每个组的观察数很少，可使用残差检查整体的正态性。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。
• 交互作用判断：两因素方差分析时，需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。当不存在交互作用时，说明两因素的作用彼此独立，逐一分析各因素的主效应即可；计算主效应时，在模型中仍需要保留交互项。