R语言学习系列27方差分析Word文档下载推荐.docx
《R语言学习系列27方差分析Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《R语言学习系列27方差分析Word文档下载推荐.docx(23页珍藏版)》请在冰豆网上搜索。
![R语言学习系列27方差分析Word文档下载推荐.docx](https://file1.bdocx.com/fileroot1/2023-1/9/01d46759-e65d-4f05-80c4-233b1d2a45a8/01d46759-e65d-4f05-80c4-233b1d2a45a81.gif)
C)。
当交互影响项呈现统计不显著时,表明各个因素独立,当呈现统计显著时,就需要列出这个交互影响项的效应,以助于作出正确的统计推断。
举例解释上述概念:
要考察焦虑症的治疗疗效,一个因素是治疗方案,有2种治疗方案,即该因素有2个水平;
(治疗方案称为组间因子,因为每个患者只能被分配到一个组别中,没有患者同时接受两种治疗);
再考虑一个因素治疗时间,也有两个水平:
治疗5周和治疗6个月,同一患者在5周和6个月不止一次地被测量(两次),称为重复测量(治疗时间称为组内因子,因为每个患者在所有水平下都进行了测量)。
建立方差分析模型时,既要考虑两个因素治疗方案和治疗时间(主效应),又要考虑治疗方案和时间的交互影响(交互效应),此时即两因素混合模型方差分析。
当某个因素的各个水平下的因变量的均值呈现统计显著性差异时,必要时可作两两水平间的比较,称为均值间的两两比较。
二、R语言实现方差分析对数据的要求:
满足正态性(来自同一正态总体)和方差齐性(各组方差相等),在这两个条件下,若各组有差异,则只可能是来自影响因素的不同水平。
用aov()函数进行方差分析,基本格式为:
aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)其中,formula为方差分析公式;
data为数据框;
projection设置是否返回预测结果;
qr设置是否返回QR分解结果;
contrasts为公式中一些因子的列表
formula公式的表示:
(y为因变量,ABC为分组因子)
符号
用法
分隔符号,左边为响应变量,右边为解释变量
eg:
y~A+B+C
+
分隔解释变量
:
表示变量的交互项eg:
y~A+B+A:
B
*
表示所有可能交互项
y~A*B*C可展开为:
y~A+B+C+A:
B+A:
C+B:
C+A:
B:
C
^
表示交互项达到次数
y~(A+B+C)^2展开为:
.
表示包含除因变量外的所有变量eg:
若一个数据框包括变量y,A、B和C,代码y~.可展开为y~A+B+C
常见研究设计的表达式:
(小写字母表示定量变量,大写字母表
示组别因子,Subject是对被试者独有的标识变量)
设计
表达式
单因素ANOVA
y~A
含单个协变量的单因素ANCOVA
y~x+A
双因素ANOVA
y~A*B
含两个协变量的双因素ANCOVA
y~x1+x2+A*B
随机化区组
y~B+A,B为区组因子
单因素组内ANOVA
y~A+Error(Subject/A)
含单个组内因子(W)和单个组间因子
(B)的重复测量ANOVA
y~B*W+Error(Subject/W)
注意:
非均衡设计时或存在协变量时,效应项的顺序对结果影响较大,越基础的效应越需要放在表达式前面,首先是协变量、然后是主效应、接着是双因素的交互项,再接着是三因素的交互项。
若研究不是正交的,一定要谨慎设置效应的顺序。
有三种类型的方法可以分解y~A+B+A:
B右边各效应对y所解释的方差:
类型I(序贯型)
效应根据表达式中先出现的效应做调整。
A不做调整,B根据A调整,A:
B交互项根据A和B调整。
类型II(分层型)
效应根据同水平或低水平的效应做调整。
A根据B调整,B依据A调整,A:
B交互项同时根据A和B调整。
类型III(边界型)
每个效应根据模型其他各效应做相应调整。
A根据B和A:
B做调整,A:
R默认调用类型I方法,其他软件(比如SAS和SPSS)默认调用类型III方法。
car包中的Anova()函数(不要与标准anova()函数混淆)提供了使用类型II或类型III方法的选项,而aov()函数使用的是类型I方法。
若想使结果与其他软件(如SAS和SPSS)提供的结果保持一致,可以使用Anova()函数。
三、单因素方差分析
1个因变量,1个影响因素:
总差异Yij=平均差异μ+因素差异αi+随机差异εij例1比较4种品牌的胶合板的耐磨性,各抽取5个样品,相同转速磨损相同时间测得磨损深度(mm),比较4个品牌胶合板的耐磨性有无差异?
部分数据如下(ex27_ex1.Rdata):
setwd("
E:
/办公资料/R语言/R语言学习系列/codes"
)load("
ex27_ex1.Rdata"
)head(datas)
wearbrand
12.30A
22.32A
32.40A
42.45A
52.58A
62.35Battach(datas)table(brand)#各组的样本数
brand
ABCD
5555
aggregate(wear,by=list(brand),mean)#各组均值
Group.1x
1A2.410
2B2.404
3C2.046
4
#各组标准差
D2.572aggregate(wear,by=list(brand),sd)
1A0.11269428
2B0.11760102
3C0.11216060
4D0.03271085
library(car)
#用Q-Q图检验数据
qqPlot(lm(wear~brand,data=datas),simulate=TRUE)的正态性
leveneTest(wear~as.factor(brand),data=datas)#方差齐性检验
Levene'
sTestforHomogeneityofVariance(center=median)
DfFvaluePr(>
F)group30.69870.5664
16
fit<
-aov(wear~brand,data=datas)#单因素方差分析(检验组间差异)
summary(fit)
DfSumSqMeanSqFvaluePr(>
F)
brand30.73980.2466024.553.15e-06***
Residuals160.16070.01005
Signif.codes:
0‘***'
0.001‘**'
0.01‘*'
0.05‘.'
0.1‘'
1
说明:
方差齐性检验,原假设H0:
方差齐,p值=0.5664>
0.05,故接受原假设,即方差齐。
单因素方差分析结果,brand是因素,Residuals是残差,各列依次为自由度、平方和、均方和、F统计量,p值=3.15e-06<
0.05,拒绝原假设,即不同品牌的磨损(均值)有显著差别。
library(gplots)
plotmeans(wear~brand,xlab="
品牌"
ylab="
磨损"
)#图形展示带95%置信区间的各组均值
通过前面的分析知道,不同品牌的磨损(均值)有显著差别,但
并不知道哪个品牌与其它品牌有显著差别。
TukeyHSD()函数提供了对
各组均值差异的成对检验
TukeyHSD(fit)
Tukeymultiplecomparisonsofmeans
95%family-wiseconfidencelevel
Fit:
aov(formula=wear~brand,data=datas)$brand
difflwruprpadj
B-A-0.006-0.187353450.17535350.9996826
C-A-0.364-0.54535345-0.18264650.0001610
D-A0.162-0.019353450.34335350.0886142
C-B-0.358-0.53935345-0.17664650.0001929
D-B0.168-0.013353450.34935350.0744337
D-C0.5260.344646550.70735350.0000019
可以看出(H0:
无差异),B与A的差异非常不显著,C与A、
C与B、D与C的差异非常显著
multcomp包中的glht()函数提供了更为全面的多重均值比较方法。
library(multcomp)attach(datas)
tuk<
-glht(fit,linfct=mcp(brand="
Tukey"
))
#注意datas$brand必须是因子型summary(tuk)
SimultaneousTestsforGeneralLinearHypotheses
MultipleComparisonsofMeans:
TukeyContrasts
0.1‘'
aov(formula=wear~brand,data=datas)
LinearHypotheses:
EstimateStd.ErrortvaluePr(>
|t|)
B-
A=
=0
-0.00600
0.06339
-0.095
0.9997
-A=
-0.36400
-5.742
<
0.001***
D
0.16200
2.556
0.0886.
-B=
-0.35800
-5.648
0.16800
2.650
0.0743.
-C=
0.52600
8.298
0.001‘**'
0.01‘*'
0.05
(Adjustedpvaluesreported--single-stepmethod)
plot(cld(tuk,level=0.05),col="
lightgrey"
)
说明:
标记相同字母(标记b)的品牌ABD认为是无显著差异,在同一亚组,而品牌C(标记a)与另外三个品牌有显著差异。
另外,也可以进行多重t检验,使用函数:
pairwise.t.test(x,g,p.adjust.method=,...)
其中,x为因变量,g为因子型的分组变量;
p.adjust.method设置p值的修正方法,由于多次重复t检验会大大增加犯第一类错误的概率,为此要进行p值的修正,使用bonferroni法修正效果较好。
pairwise.t.test(wear,brand,p.adjust.method="
bonferroni"
PairwisecomparisonsusingttestswithpooledSDdata:
wearandbrand
ABC
B1.00000--
C0.000180.00022-
D0.126950.104742.1e-06
Pvalueadjustmentmethod:
bonferroni
原假设H0:
无差异,可见A与B无差异,C与ABD有显著差异。
最后,方差分析对离群点非常敏感,检验是否有离群点:
outlierTest(fit)
NoStudentizedresidualswithBonferonnip<
Largest|rstudent|:
rstudentunadjustedp-valueBonferonnip
92.5281030.0231820.46364
经检验无离群点。
三、两因素方差分析
1个因变量,2个影响因素:
总差异Yijk=平均差异μ+因素1差异αi+因素2差异βi
+因素1,2交互作用差异γij+随机差异εijk
例2研究60只豚鼠的牙齿生长数据,按2种喂食方法:
橙汁、维生素C,各喂食方法中抗坏血酸含量都有3个水平:
0.5mg、1mg、2mg,
分配为6组,每组各10只,牙齿长度为因变量。
做两因素方差分析。
attach(ToothGrowth)head(ToothGrowth)
lensuppdose
14.2VC0.5
211.5VC0.5
37.3VC0.5
45.8VC0.5
56.4VC0.5
610.0VC0.5
table(supp,dose)#各组样本数相同,即为均衡设计
dose
supp0.512
OJ101010
VC101010aggregate(len,by=list(supp,dose),mean)
Group.1Group.2x
#计算各组均值
1
OJ
0.513.23
2
VC
0.57.98
3
1.022.70
1.016.77
5
2.026.06
6
2.026.14
#计算各组标准差
0.54.459709
0.52.746634
1.03.910953
1.02.515309
2.02.655058
2.04.797731
aggregate(len,by=list(supp,dose),sd)
bartlett.test(len~supp,data=ToothGrowth)检验
#关于因素supp的方差齐性
Bartletttestofhomogeneityofvariances
data:
lenbysupp
Bartlett'
sK-squared=1.4217,df=1,p-value=0.2331bartlett.test(len~dose,data=ToothGrowth)#关于因素dose的方差齐性
检验
lenbydose
sK-squared=0.66547,df=2,p-value=0.717fit<
-aov(len~supp*dose,data=ToothGrowth)#做两因素方差分析,考虑全部效应
supp1205.4205.412.3170.000894***
dose12224.32224.3133.415<
2e-16***
Signif.codes:
可以看出,主效应supp和dose都非常显著(p值都远小
于0.05),交互效应也显著(p值=0.0246<
0.05)。
若交互作用不显著,
可以可以只做去掉交互效应的方差分析。
图形化展示两因素方差分析的交互效应:
par(mfrow=c(1,2))
interaction.plot(dose,supp,len,type="
b"
col=c("
red"
"
blue"
),pch=c(16,18),main="
InteractionbetweenDoseandSupp"
interaction.plot(supp,dose,len,type="
col=c("
InteractionbetweenSuppandDose"
有一个图的线有交叉,说明有交互作用。
可以看出随着橙汁和维生素C中的抗坏血酸剂量的增加,牙齿长度变长;
。
对于0.5mg和1mg剂量,橙汁比维生素C更能促进牙齿生长;
对于2mg剂量的
抗坏血酸,两种喂食方法下牙齿长度增长相同
也可以用HH包中的interaction2wt()函数(也适合三因素方差分
析)来展示更全面的可视化结果:
library(HH)
interaction2wt(len~supp*dose)
三、重复测量方差分析
重复测量方差分析,即受试者被测量不止一次
例3(1个组内1个组间因子的重复测量)在某浓度CO2的环境中,
对寒带植物(来自魁北克)和非寒带植物的(来自密西西比)光合作
用率进行比较。
因变量uptake为CO2吸收量,自变量Type(组间因
子)为植物类型,自变量conc(组内因子)为七种水平的CO2浓度。
attach(CO2)
head(CO2)#注意CO2是长格式的数据
PlantTypeTreatmentconcuptake
1Qn1Quebecnonchilled9516.0
2Qn1Quebecnonchilled17530.4
3Qn1Quebecnonchilled25034.8
4Qn1Quebecnonchilled35037.2
5Qn1Quebecnonchilled50035.3
6Qn1Quebecnonchilled67539.2w1b1<
-subset(CO2,Treatment=="
chilled"
)#先只考虑寒带植物
-aov(uptake~(conc*Type)+Error(Plant/conc),data=w1b1)summary(fit)
Error:
Plant
Type12667.22667.260.410.00148**
Residuals4176.644.1
Plant:
conc
conc1888.6888.6215.460.000125***
conc:
Type1
Residuals4
239.2
16.5
4.1
58.010.001595**
0.001‘**'
0.01
‘*'
0.