R语言学习系列27方差分析.docx

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

---

: