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

发布于 2021年12月29日 星期三 13:10:00 浏览:6566
原创不易,转载请注明来源,感谢!
附件下载:
两因素方差分析二.zip 请勿重复点击,如无响应请耐心等待或稍后再试。

在前面文章中介绍了两因素方差分析(Two-way ANOVA)——不存在交互作用时在SAS软件中的实现,本篇文章将实例演示在SAS软件中实现两因素方差分析——存在交互作用时的操作步骤。

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

一、案例介绍

观察A、B两种镇痛药物联合运用在产妇分娩时的镇痛效果。A药(Drug_A)取3个剂量:1.0、2.5、5.0 mg;B药(Drug_B)也取3个剂量:5.0、15.0、30.0μg,共9个处理组。将27名产妇随机分成9组,每组3名产妇,记录每名产妇分娩时的镇痛时间Time (min)。试分析A、B两药联合运用的镇痛效果。

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

图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\两因素方差分析2.csv "
dbms=csv replace;
getnames=yes;
run;

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

(二) 适用条件判断

1. 生成因变量残差

(1) 软件操作

①运用SORT过程以Drug_A、Drug_B升序排列的形式重新排列,具体代码如下所示:

Proc sort data=mydata.example;
by Drug_A Drug_B;
run;

重新排列后数据如图2所示:

图2

②运用REG过程计算残差,具体代码如下所示:

Proc reg data=mydata.example;
by Drug_A Drug_B;
model time= Drug_A Drug_B /r;
output out=mydata.example2 r=residual; 
run; 
quit;

BY语句以“Drug_A”和“Drug_B”变量分组计算残差,MODEL语句以“time”为因变量,“Drug_A”和“Drug_B”作为分组变量计算残差的值。OUTPUT语句将残差输入新建的表“mydata.example2”,新建变量命名为“residual”。“r=residual”将该步骤中计算得的残差输入至新表。

(2) 结果解读

上述操作将可以得到生成的新变量“residual”,如图3所示。

图3

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

(1) 软件操作

①运用UNIVARIATE语句进行检验,具体代码如下所示:

proc univariate data=mydata.example2 normal plot;
class Drug_A Drug_B;
var residual;
run;

CLASS语句起到拆分文件的作用,分别分析“Drug_A”和“Drug_B”分别为“1”“2”“3”时共9个处理组是否存在异常值。PLOT选项可作各分组的平行条状图,箱线图和正态概率图。“Drug_A”=1且“Drug_B”=1时如图4所示,其余图像不一一列举。

图4

②运用ANOVA过程作各组箱式图,具体代码如下所示:

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

该步骤中,CLASS语句指定分类变量,即“Drug_A”和“Drug_B”两变量。MODEL语句指定分析因素与响应变量,即按照“Drug_A”和“Drug_B”同时分组对“residual”变量分析。该过程输出结果如图5所示,图中显示的横坐标“Drug_A*Drug_B”即表示“Drug_A”和“Drug_B”分别为“1”“2”“3”时共9个处理组的箱式图。

图5
(2) 结果解读

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

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

两因素方差分析时,可分别考察每一组原始数据的正态性或使用残差考察整体正态性。

(1) 软件操作

运用UNIVARIATE语句对残差考察整体正态性,具体代码如下所示:

proc univariate data=mydata.example2 normal plot;
var residual;
run;

NORMAL选项可以对进行整体残差的正态性检验。PLOT选项可以得到平行条状图,箱线图和正态概率图。

(2) 结果解读

图6显示了四种正态性检验的结果,Shapiro-Wilk (夏皮罗-威尔克正态性,S-W)检验和Kolmogorov-Smirnov (柯尔莫哥洛夫-斯米诺夫,K-S)检验。K-S检验适用于大样本资料,本案查看S-W检验结果,可见P值>0.1,提示整体残差服从正态分布。

图6

图7为整体残差正态性检验的正态概率图,可见散点基本围绕对角线分布,也提示数据服从正态分布。综上,本案例满足条件5。

图7

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

(1) 软件操作

运用ANOVA过程进行检验,具体代码如下所示:

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

该步骤中,CLASS语句指定分类变量,即“Drug_A”和“Drug_B”两变量。MODEL语句指定分析因素与响应变量,即按照“Drug_A”和“Drug_B”同时分组对“time”变量分析。 MEANS语句计算均数和方差,其常见选项HOVTEST可指定分析方法,此处我们使用LEVENE检验,“(type=abs)”指定该过程计算绝对值(absolute value)进行分析。

(2) 结果解读

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

图8

(三) 统计描述及推断

1. 软件操作

①前一步骤中UNIVARIATE语句可得各组数据的统计描述结果。也可用MEANS过程,其代码如下:

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

②运用ANOVA过程进行检验,该过程与上述ANOVA过程稍有不同,具体代码如下所示:

proc anova data = mydata.example2;
class Drug_A Drug_B;
model time = Drug_A|Drug_B;
means Drug_A|Drug_B /snk;
run;

该过程中使用“Drug_A|Drug_B”,可以对两变量的交互作用进行检验,其也可以替换为“Drug_A Drug_B Drug_A*Drug_B”。SNK法指定两两比较方法。

2. 结果解读

(1) 统计描述

图9是由MEANS过程所得的“Descriptives Statistics(描述统计)”列出了各组的均值和标准差,可知A药在1.0 mg水平,B药为5.0、15.0、30.0 μg水平时,镇痛时间分别为:83.33±20.21、100.00±18.03、85.00±10.00 min;A药在2.5 mg水平,B药为5、15、30 μg水平时,镇痛时间分别为:90.00±21.79、115.00±21.79、135.00±15.00 min;A药在5.0 mg水平,B药为5、15、30 μg水平时,镇痛时间分别为:110.00±21.79、95.00±27.84、176.67±15.28 min。ANOVA过程给出time分布的箱线图,如图10所示。

图9
图10
(2) 统计推断

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

图11
(3) 交互作用判断

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

图10第三行是对A药和B药交互作用的判断,FDrug_A*Drug_B=5.07,P值为0.0065,有统计学意义,说明两者存在交互作用。

(四) 单独效应分析

1. 软件操作

proc glm data = mydata.example2;
class Drug_A Drug_B;
model time = Drug_A Drug_B Drug_A*Drug_B;
lsmeans Drug_A *Drug_B /slice=Drug_B;
run;
quit;

更改lsmeans Drug_A /slice=Drug_B中变量名顺序,对不同的A药浓度下分析B药的作用,具体代码如下:

proc glm data = mydata.example2;
class Drug_A Drug_B;
model time = Drug_A Drug_B Drug_A*Drug_B;
lsmeans Drug_A*Drug_B /slice=Drug_A;
run;
quit;

2. 结果解读

(1) A药的单独效应

A药的单独效应就是在不同的B药浓度下分析A药的作用。由图10可大致看出,在不同的B药浓度下,A药对镇痛时间的影响不同,具体的检验结果见图12。结果显示,在药物B为5.0 μg水平时,药物A各水平之间的镇镇痛时间差异无统计学意义(F=1.49,P=0.252);在药物B为15.0 μg水平时,药物A各水平之间的镇镇痛时间差异无统计学意义(F=0.84,P=0.449);但在药物B为30.0 μg水平时,药物A各水平之间的镇镇痛时间差异有统计学意义(F=16.29,P<0.001)。

图12
(2) B药的单独效应

由图10可大致看出,在不同的A药浓度下,B药对镇痛时间的影响不同,具体的检验结果见图13。结果显示,在药物A为1.0 mg水平时,药物B各水平之间的镇痛时间差异无统计学意义(F=0.65,P=0.533);在药物A为2.5 mg水平时,药物B各水平之间的镇痛时间差异有统计学意义(F=3.93,P=0.038);但在药物A为5.0 mg水平时,药物B各水平之间的镇痛时间差异有统计学意义(F=14.61,P<0.001)。

图13

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

1. 软件操作

ods output Estimates=estimate;
proc glm data = mydata.example2;
    class Drug_A Drug_B;
    model time =Drug_A Drug_B Drug_A*Drug_B;
	 estimate '(Drug_A 1.0mg vs 2.5mg)/Drug_B=30ug'    Drug_A 1 -1 0  Drug_A*Drug_B 0 0 1 0 0 -1 0 0 0;
	 estimate '(Drug_A 1.0mg vs 5.0mg)/Drug_B=30ug'    Drug_A 1 0 -1  Drug_A*Drug_B 0 0 1 0 0 0 0 0 -1;
	 estimate '(Drug_A 2.5mg vs 5.0mg)/Drug_B=30ug'    Drug_A 0 1 -1  Drug_A*Drug_B 0 0 0 0 0 1 0 0 -1;
	 estimate '(Drug_B 5.0ug vs 15.0ug)/Drug_A=2.5mg'  Drug_B 1 -1 0  Drug_A*Drug_B 0 0 0 1 -1 0 0 0 0;
	 estimate '(Drug_B 5.0ug vs 30.0ug)/Drug_A=2.5mg'  Drug_B 1 0 -1  Drug_A*Drug_B 0 0 0 1 0 -1 0 0 0;
	 estimate '(Drug_B 15.0ug vs 30.0ug)/Drug_A=2.5mg' Drug_B 0 1 -1  Drug_A*Drug_B 0 0 0 0 1 -1 0 0 0;
	 estimate '(Drug_B 5.0ug vs 15.0ug)/Drug_A=5mg'    Drug_B 1 -1 0  Drug_A*Drug_B 0 0 0 0 0 0 1 -1 0;
	 estimate '(Drug_B 5.0ug vs 30.0ug)/Drug_A=5mg'    Drug_B 1 0 -1 Drug_A*Drug_B 0 0 0 0 0 0 1 0 -1;
	 estimate '(Drug_B 15.0ug vs 30.0ug)/Drug_A=5mg'   Drug_B 0 1 -1  Drug_A*Drug_B 0 0 0 0 0 0 0 1 -1; 
run;
quit;
proc print data=estimate;run;
data estimate;
set estimate;
Probt=Probt*3;/* 通过“Bonferroni (邦弗罗尼)”法两两比较*/
proc print;run;

2. 结果解读

图14

图14为B药30 µg浓度下,A药的两两检验结果和A药在2.5 mg和5 mg情况下,B药的两两检验结果,提供了不同浓度间两两比较的“LSMEAN (最小二乘均值)”和各组相比的“Pr (P值)”。

由之前结果可知,当B药为5 μg与15 μg 时,A药整体无统计学意义(P>0.05),故这里不再分析。当B药为30 μg时, A药浓度为2.5 mg和5.0 mg与1.0 mg相比,均值差逐渐增大。2.5 mg比1.0 mg时均值增加了50.00 min,差异有统计学意义(P=0.0182);至5.0 mg时增加了91.67 min,差异有统计学意义(P<0.0001);2.5 mg到5.0 mg时增加了41.67 min,差异无统计学意义(P=0.0554)。表明,在B药为30.0 μg水平时,随着A药浓度的增加,镇痛时间呈上升趋势。

当A药为1.0 mg时,B药整体无统计学意义(P>0.05),这里不再比较。当药物A为2.5 mg时,药物B浓度30.0 μg与5.0 μg相比,镇痛时间差异有统计学意义(P=0.036),30.0 μg与15.0 μg相比以及15.0 μg与5.0 μg相比差异均无统计学意义(P>0.05)。当药物A为5.0 mg时,药物B浓度30.0 μg与15.0 μg相比以及30.0 μg与5.0 μg相比,镇痛时间差异均有统计学意义(P<0.01),但15.0 μg与5.0 μg相比,镇痛时间差异无统计学意义(P=1.090)。

四、结论

本研究采用两因素方差分析探讨A、B两种镇痛药物联合运用在产妇分娩时的镇痛效果。通过对模型残差绘制箱线图提示,数据不存在异常值;通过Shapiro-Wilk检验,提示各组数据服从正态分布;通过Levene’s检验,提示数据总体方差相等;两药物之间存在交互作用(FDrug_A*Drug_B=5.073,P=0.006),故进行简单效应分析。

药物A在1.0 mg水平,药物B为5.0、15.0、30.0 μg水平时,镇痛时间分别为:83.333±20.207、100.000±18.028、85.000±10.000 min;药物A在2.5 mg水平,药物B为5、15、30 μg水平时,镇痛时间分别为:90.000±21.794、115.000±21.794、135.000±15.000 min;药物A在5.0 mg水平,药物B为5、15、30 μg水平时,镇痛时间分别为:110.000±21.794、95.000±27.839、176.667±15.275 min。

药物A的简单效应分析显示,在药物B为30.0 μg水平时,药物A各水平之间的镇痛时间差异有统计学意义(F=16.29,P<0.001)。但在药物B为5.0 μg和15.0 μg水平时,药物A各水平之间的镇痛时间差异均无统计学意义(P>0.05)。对药物A在药物B为30.0 μg水平时进行“Bonferroni (邦弗罗尼)”法两两比较结果显示, 1.0 mg至2.5 mg时均值增加了50.00 min,至5.0 mg时均值增加了91.67 min,差异均有统计学意义(P<0.001);表明,在药物B为30.0 μg水平时,随着药物A浓度的增加,镇痛时间呈上升趋势。

药物B的简单效应分析显示,在药物A为1.0 mg水平时,药物B各水平之间的镇痛时间差异无统计学意义(F=0.65,P=0.533);但在药物A为2.5 mg和5.0 mg水平时,药物B各水平之间的镇痛时间差异均有统计学意义(P<0.05)。“Bonferroni (邦弗罗尼)”法两两比较显示,当药物A为2.5 mg时,药物B浓度30.0 μg与5.0 μg相比,镇痛时间差异有统计学意义(P=0.036),30.0 μg与15.0 μg相比以及15.0 μg与5.0 μg相比差异均无统计学意义(P>0.05)。当药物A为5.0 mg时,药物B浓度30.0 μg与15.0 μg相比以及30.0 μg与5.0 μg相比,镇痛时间差异均有统计学意义(P<0.01),但15.0 μg与5.0 μg相比,镇痛时间差异无统计学意义(P=1.000)。

综上可知,对产妇分娩时镇痛,A、B两种药物联合运用时会相互影响,随着药物浓度的增加,联合作用趋于复杂和不稳定。

五、分析小技巧

  • 交互作用判断:两因素方差分析时,需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义,则需要分析单独效应。此时,单纯研究某个因素的作用并无意义,应分别探讨另一个因素不同水平时对该因素的作用。当不存在交互作用时,说明两因素的作用彼此独立,逐一分析各因素的主效应即可;计算主效应时,在模型中仍需要保留交互项。
End
文章目录 沉浸式阅读