R语言学习系列27方差分析.docx
《R语言学习系列27方差分析.docx》由会员分享,可在线阅读,更多相关《R语言学习系列27方差分析.docx(25页珍藏版)》请在冰豆网上搜索。
R语言学习系列27方差分析
R语言学习系列27-方差分析
22.方差分析
一、方差分析原理
1.方差分析概述
方差分析可用来研究多个分组的均值有无差异,其中分组是按影响因素的不同水平值组合进行划分的。
方差分析是对总变异进行分析。
看总变异是由哪些部分组成的,这些部分间的关系如何。
方差分析,是用来检验两个或两个以上均值间差别显著性(影响观察结果的因素:
原因变量(列变量)的个数大于2,或分组变量(行变量)的个数大于1)。
一元时常用F检验(也称一元方差分析),多元时用多元方差分析(最常用Wilks’∧检验)。
方差分析可用于:
(1)完全随机设计(单因素)、随机区组设计(双因素)、析因设计、拉丁方设计和正交设计等资料;
(2)可对两因素间交互作用差异进行显著性检验;
(3)进行方差齐性检验。
要比较几组均值时,理论上抽得的几个样本,都假定来自正态总体,且有一个相同的方差,仅仅均值可以不相同。
还需假定每一个观察值都由若干部分累加而成,也即总的效果可分成若干部分,而每一部分都有一个特定的含义,称之谓效应的可加性。
所谓的方差是离均差平方和除以自由度,在方差分析中常简称为均方(MeanSquare)。
2.基本思想
基本思想是,将所有测量值上的总变异按照其变异的来源分解为多个部份,然后进行比较,评价由某种因素所引起的变异是否具有统计学意义。
根据效应的可加性,将总的离均差平方和分解成若干部分,每一部分都与某一种效应相对应,总自由度也被分成相应的各个部分,各部分的离均差平方除以各自的自由度得出各部分的均方,然后列出方差分析表算出F检验值,作出统计推断。
方差分析的关键是总离均差平方和的分解,分解越细致,各部分的含义就越明确,对各种效应的作用就越了解,统计推断就越准确。
效应项与试验设计或统计分析的目的有关,一般有:
主效应(包括各种因素),交互影响项(因素间的多级交互影响),协变量(来自回归的变异项),等等。
当分析和确定了各个效应项S后,根据原始观察资料可计算出各个离均差平方和SS,再根据相应的自由度df,由公式MS=SS/df,求出均方MS,最后由相应的均方,求出各个变异项的F值,F值实际上是两个均方之比值,通常情况下,分母的均方是误差项的均方。
根据F值的分子、分母均方的自由度f1和f2,在确定显著性水平为α情况下,由F(f1,f2)临界值表查得单侧Fα界限值。
当Fα,不拒绝原假设H0,说明不拒绝这个效应项的效应为0的原假设,也即这个效应项是可能对总变异没有实质影响的;若F>Fα则P值≤α,拒绝原假设H0,也即这个效应项是很可能对总变异有实质影响的。
3.方差分析的实验设计
为了确定方差分析表中各个有关效应项,需要在试验设计阶段就作出安排,再根据设计要求进行试验,得出原始观察值,按原来设计方案算出方差分析表中的各项。
在试验设计阶段通常需要考虑如下4个方面:
(1)研究的因变量
即试验所要观察的主要指标,一次试验时可以有多个观察指标,方差分析时也可以同时对多个因变量进行分析;
(2)因素和水平
试验的因素(factor)可以是品种、人员、方法、时间、地区等等,因素所处的状态叫水平(level)。
在每一个因素下面可以分成若干水平。
(3)因素间的交互影响
多因素的试验设计,有时需要分析因素间的交互影响(interaction),2个因素间的交互影响称为一级交互影响(A×B);3个因素间的交互影响称为二级交互影响(A×B×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
*
表示所有可能交互项
eg:
y~A*B*C可展开为:
y~A+B+C+A:
B+A:
C+B:
C+A:
B:
C
^
表示交互项达到次数
eg:
y~(A+B+C)^2展开为:
y~A+B+C+A:
B+A:
C+B:
C
.
表示包含除因变量外的所有变量
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:
B交互项根据A和B调整。
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个品牌胶合板的耐磨性有无差异部分数据如下():
setwd("E:
/办公资料/R语言/R语言学习系列/codes")
load("")
head(datas)
wearbrand
1A
2A
3A
4A
5A
6B
attach(datas)
table(brand)#各组的样本数
brand
ABCD
5555
aggregate(wear,by=list(brand),mean)#各组均值
x
1A
2B
3C
4D
aggregate(wear,by=list(brand),sd)#各组标准差
x
1A0.
2B0.
3C0.
4D
library(car)
qqPlot(lm(wear~brand,data=datas),simulate=TRUE)#用Q-Q图检验数据的正态性
leveneTest(wear~(brand),data=datas)#方差齐性检验
Levene'sTestforHomogeneityofVariance(center=median)
DfFvaluePr(>F)
group3
16
fit<-aov(wear~brand,data=datas)#单因素方差分析(检验组间差异)
summary(fit)
DfSumSqMeanSqFvaluePr(>F)
brand3***
Residuals16
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
说明:
方差齐性检验,原假设H0:
方差齐,p值=>,故接受原假设,即方差齐。
单因素方差分析结果,brand是因素,Residuals是残差,各列依次为自由度、平方和、均方和、F统计量,p值=<,拒绝原假设,即不同品牌的磨损(均值)有显著差别。
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.
C-A-0.
D-A
C-B-0.
D-B
D-C0.
说明:
可以看出(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
Fit:
aov(formula=wear~brand,data=datas)
LinearHypotheses:
EstimateStd.ErrortvaluePr(>|t|)
B-A==0
C-A==0<***
D-A==0.
C-B==0<***
D-B==0.
D-C==0<***
---
Signif.codes:
0‘***’‘**’‘*’‘.’‘’1
(Adjustedpvaluesreported--single-stepmethod)
plot(cld(tuk,level=,col="lightgrey")
说明:
标记相同字母(标记b)的品牌ABD认为是无显著差异,在同一亚组,而品牌C(标记a)与另外三个品牌有显著差异。
另外,也可以进行多重t检验,使用函数:
g,其中,x为因变量,g为因子型的分组变量;
设置p值的修正方法,由于多次重复t检验会大大增加犯第一类错误的概率,为此要进行p值的修正,使用bonferroni法修正效果较好。
"bonferroni")
PairwisecomparisonsusingttestswithpooledSD
data:
wearandbrand
ABC
B--
C-
D
Pvalueadjustmentmethod:
bonferroni
说明:
原假设H0:
无差异,可见A与B无差异,C与ABD有显著差异。
最后,方差分析对离群点非常敏感,检验是否有离群点:
library(car)
outlierTest(fit)
NoStudentizedresidualswithBonferonnip<
Largest|rstudent|:
rstudentunadjustedp-valueBonferonnip
9
说明:
经检验无离群点。
三、两因素方差分析
1个因变量,2个影响因素:
总差异Yijk=平均差异μ+因素1差异αi+因素2差异βi
+因素1,2交互作用差异γij+随机差异εijk
例2研究60只豚鼠的牙齿生长数据,按2种喂食方法:
橙汁、维生素C,各喂食方法中抗坏血酸含量都有3个水平:
、1mg、2mg,分配为6组,每组各10只,牙齿长度为因变量。
做两因素方差分析。
attach(ToothGrowth)
head(ToothGrowth)
lensuppdose
1VC
2VC
3VC
4VC
5VC
6VC
table(supp,dose)#各组样本数相同,即为均衡设计
dose
supp12
OJ101010
VC101010
aggregate(len,by=list(supp,dose),mean)#计算各组均值
x
1OJ
2VC
3OJ
4VC
5OJ
6VC
aggregate(len,by=list(supp,dose),sd)#计算各组标准差
x
1OJ
2VC
3OJ
4VC
5OJ
6VC
(len~supp,data=ToothGrowth)#关于因素supp的方差齐性检验
Bartletttestofhomogeneityofvariances
data:
lenbysupp
Bartlett'sK-squared=,df=1,p-value=
(len~dose,data=ToothGrowth)#关于因素dose的方差齐性检验
Bartletttestofhomogeneityofvariances
data:
lenbydose
Bartlett'sK-squared=,df=2,p-value=
fit<-aov(len~supp*dose,data=ToothGrowth)
#做两因素方差分析,考虑全部效应
summary(fit)
DfSumSqMeanSqFvaluePr(>F)
supp1***
dose1<2e-16***
supp:
dose1*
Residuals56
---
:
0‘***’‘**’‘*’‘.’‘’1
说明:
可以看出,主效应supp和dose都非常显著(p值都远小于),交互效应也显著(p值=<)。
若交互作用不显著,可以可以只做去掉交互效应的方差分析。
图形化展示两因素方差分析的交互效应:
par(mfrow=c(1,2))
(dose,supp,len,type="b",col=c("red","blue"),pch=c(16,18),main="InteractionbetweenDoseandSupp")
(supp,dose,len,type="b",col=c("red","blue"),pch=c(16,18),main="InteractionbetweenSuppandDose")
说明:
有一个图的线有交叉,说明有交互作用。
可以看出随着橙汁和维生素C中的抗坏血酸剂量的增加,牙齿长度变长;。
对于mg
和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
1Qn1Quebecnonchilled95
2Qn1Quebecnonchilled175
3Qn1Quebecnonchilled250
4Qn1Quebecnonchilled350
5Qn1Quebecnonchilled500
6Qn1Quebecnonchilled675
w1b1<-subset(CO2,Treatment=="chilled")#先只考虑寒带植物
fit<-aov(uptake~(conc*Type)+Error(Plant/conc),data=w1b1)
summary(fit)
Error:
Plant
DfSumSqMeanSqFvaluePr(>F)
Type1**
Residuals4
---
:
0‘***’‘**’‘*’‘.’‘’1
Error:
Plant:
conc
DfSumSqMeanSqFvaluePr(>F)
conc1***
conc:
Type1**
Residuals4
---
:
0‘***’‘**’‘*’‘.’‘’1
Error:
Within
DfSumSqMeanSqFvaluePr(>F)
Residuals30869
说明:
在的显著水平下,主效应“类型”(p值=)和“浓度”(p值=)以及交叉效应“类型*浓度”(p值=)都非常显著。
attach(w1b1)
(conc,Type,uptake,type="b",col=c("red","blue"),pch=c(16,18),main="InteractionPlotforPlantTypeandConcentration")
boxplot(uptake~Type*conc,data=w1b1,col=(c("gold","green")),main="ChilledQuebecandMississippiPlants",ylab="Carbondioxideuptakerate(umol/m^2sec)")
detach(w1b1)
可以看出,魁北克省的植物比密西西比州的植物二氧化碳吸收率高,而且随着CO2浓度的升高,差异越来越明显。
注1:
重复测量设计时,需要有长格式数据才能拟合模型,若是宽数据,需要用reshape2包中的melt()函数转化为长数据;
注2:
这里是用的传统的重复测量方差分析,假设任意组内因子的协方差矩阵为球形,并且任意组内因子两水平间的方差之差都相等。
实际中,该假设一般不满足,可以尝试:
使用lme4包中的lmer()函数拟合线性混合模型(Bates,2005);
使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足(例如Geisser-Greenhouse校正);
使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型(UCLA,2009);
用多元方差分析对重复测量数据进行建模(Hand,1987)。
四、多元方差分析
1.因变量不只一个的方差分析,称为多元方差分析。
要求数据满足:
多元正态性、方差—协方差矩阵同质性(即指各组的协方差矩阵相同,通常可用Box’sM检验来评估该假设)。
例4使用MASS包中的UScereal数据集,研究美国谷物中的卡路里、脂肪、糖含量是否会因为储物架位置的不同而发生变化。
自变量货架位置shelf有1,2,3三个水平。
library(MASS)
attach(UScereal)
y<-cbind(calories,fat,sugars)
head(y)
caloriesfatsugars
[1,]
[2,]
[3,]
[4,]
[5,]
[6,]
aggregate(y,by=list(shelf),mean)
caloriesfatsugars
11
22
33
cov(y)#求协方差矩阵,主对角线为方差
caloriesfatsugars
calories
fat
sugars
library(mvnormtest)
(t(y))#多元正态性检验,注意要对y转置
Shapiro-Wilknormalitytest
data:
Z
W=,p-value=
fit<-manova(y~shelf)#做多元方差分析
summary(fit)
DfPillaiapproxFnumDfdenDfPr(>F)
shelf1361**
Residuals63
---
:
0‘***’‘**’‘*’‘.’‘’1
(fit)#输出每个单变量的方差分析结果
Responsecalories:
DfSumSqMeanSqFvaluePr(>F)
shelf14531345313***
Residuals632039823238
---
:
0‘***’‘**’‘*’‘.’‘’1
Responsefat:
DfSumSqMeanSqFvaluePr(>F)
shelf1**
Residuals63
---
: