《统计学-基于R》(08)第8章-方差分析(R3)

数据分析

(方法与案例)

作者贾俊平版权所有违者必究StatisticswithR统计学R语言第8章方差分析8.1

方差分析的原理8.2单因子方差分析8.3双因子方差分析8.4方差分析的假定及其检验8.5单因子方差分析的非参数方法NOVAA8.1方差分析的原理

8.1.1什么是方差分析

8.1.2误差分解第8章方差分析8.1.1什么是方差分析8.1方差分析的原理2018-9-25什么是方差分析(ANOVA)

(analysisofvariance)

方差分析的基本原理是在20世纪20年代由英国统计学家RonaldA.Fisher在进行实验设计时为解释实验数据而首先引入的分析各分类自变量对数值因变量影响的一种统计方法研究分类型自变量对数值型因变量的影响

一个或多个分类型自变量两个或多个(k个)处理水平或分类一个数值型因变量有单因子方差分析和双因子方差分析单因子方差分析：涉及一个分类的自变量双因子方差分析：涉及两个分类的自变量2018-9-25什么是方差分析

(例题分析)品种1品种2品种3817176827279797277816676787278897789928187877784857387867987【例8-1】（数据：example8_1.RData）为分析小麦品种对产量的影响，一家研究机构挑选了3个小麦品种：品种1、品种2、品种3，然后选择条件和面积相同的30个地块，每个品种在10个地块上试种，实验获得的产量数据如表8-1所示。2018-9-25什么是方差分析

(例题分析)分析“小麦品种”对“产量”的影响如果只分析品种一个因子对产量的影响，则称为单因子方差分析(one-wayanalysisofvariance)如果两个因子对产量的单独影响，但不考虑它们对产量的交互效应(interaction)，则称为只考虑主效应(maineffect)的双因子方差分析，或称为无重复双因子分析(two-factorwithoutreplication)如果除了考虑两个因子对产量的单独影响外，还考虑二者对产量的交互效应，则称为考虑交互效应的双因子方差分析，或称为可重复双因子分析(two-factorwithreplication)8.1.2误差分解8.1方差分析的原理2018-9-25方差分析的基本原理

(误差分解)总误差(totalerror)反映全部观测数据的误差所抽取的全部30个地块的产量之间差异处理误差(treatmenterror)—组间误差(between-grouperror)由于不同处理造成的误差，它反映了处理(品种)对观测数据(产量)的影响，因此称为处理效应(treatmenteffect)随机误差(randomerror)—组内误差(within-grouperror)由于随机因子造成的误差，也简称为误差(error)2018-9-25方差分析的基本原理

(误差分解)数据的误差用平方和(sumofsquares)表示，记为SS总平方和(sumofsquaresfortotal)，记为SST反映全部数据总误差大小的平方和抽取的全部30个地块产量之间的误差平方和处理平方和(treatmentsumofsquares)，记为SSA反映处理误差大小的平方和也称为组间平方和(between-groupsumofsquares)误差平方和(sumofsquaresoferror)，记为SSE反映随机误差大小的平方和称为误差平方和也称为组内平方和(within-groupsumofsquares)2018-9-25方差分析的基本原理

(误差分解)2018-9-25方差分析的基本原理

(误差分析)方差分析的基本原理就是要分析数据的总误差中有没有处理误差。如果处理(超市的不同位置)对观测数据(销售额)没有显著影响，意味着没有处理误差。这时，每种处理所对应的总体均值(

i)应该相等如果存在处理误差，每种处理所对应的总体均值(

i)至少有一对不相等就例8—1而言，在只考虑品种一个因子的情况下，方差分析也就是要检验下面的假设H0：

1

2

3

H1：

1,

2,

3

不全相等8.2单因子方差分析

8.2.1数学模型

8.2.2效应检验

8.2.3效应量分析

8.2.4多重比较第8章方差分析8.2.1数学模型8.2单因子方差分析2018-9-25单因子方差分析

(数学模型)设因子A有I种处理(比如品种有“品种1”、“品种2”、“品种3”3种处理)，单因子方差分析可用下面的线性模型来表示

2018-9-25单因子方差分析

(数学模型)设全部观测数据的总均值为

，第i个处理效应用第i个处理均值与总均值的差(

i-)

表示，记为

i，即i=i-

。这样，第i个处理均值被分解成i=i+，方差分析模型可以表达为

8.2.2效应检验8.2单因子方差分析2018-9-25提出假设一般提法H0

：

i

=0

(i=1,2,…,I)没有处理效应

H1：

i

至少有一个不等于0有处理效应

注意：拒绝原假设，只表明至少有一个处理的效应显著，并不意味着所有的粗粒的效应都显著2018-9-25构造检验的统计量F单因子方差分析的方差分析表误差来源平方和SS自由度df均方MS检验统计量F处理效应误差

总效应

2018-9-25单因子方差分析

(例题分析)

#将表8-1的短格式数据转为长格式数据，并显示前6列load("C:/example/ch8/example8_1.RData")example8_1<-cbind(example8_1,id=factor(1:10))library(reshape)example8_2<-melt(example8_1,id.vars=c("id"),variable_name="品种")example8_2<-rename(example8_2,c(id="地块",value="产量"))save(example8_2,file="C:/example/ch8/example8_2.RData")example8_22018-9-25#绘制3个品种数据产量的箱线图

#计算描述统计量3个品种产量的描述性分析load("C:/example/ch8/example8_2.RData")attach(example8_2)boxplot(产量~品种,data=example8_2,col="gold",main="",ylab="产量",xlab="品种")my_summary<-function(x){with(x,data.frame("均值"=mean(产量),"标准差"=sd(产量),n=length(产量)))}library(plyr)ddply(example8_2,.(品种),my_summary)2018-9-25#方差分析表#方差分析模型的参数估计#绘制均值图方差分析表和参数估计attach(example8_2)model_1w<-aov(产量~品种)summary(model_1w)model_1w$coefficientslibrary(gplots)plotmeans(产量~品种,data=example8_2)8.2.3效应量分析8.2单因子方差分析2018-9-25效应量分析

2018-9-25library(DescTools)

model_1w<-aov(产量~品种)

EtaSq(model_1w,anova=T)效应量分析

【例8-2】的效应量8.2.3多重比较8.2单因子方差分析2018-9-25多重比较的意义通过对总体均值之间的配对比较来进一步检验到底哪些均值之间存在差异比较方法有多种，如Fisher的LSD方法、Tukey-Kramer的HSD方法等2018-9-25Fisher的LSD方法LSD是最小显著差异(leastsignificantdifference)的缩写，该检验方法是由统计学家Fisher提出来的，因此也称为Fisher的最小显著差异方法，简称LSD方法LSD的适用场合：如果研究者在事先就已经计划好要对某对或某几对均值进行比较，不管方差分析的结果如何(拒绝或不拒绝原假设)，都要进行比较，这时适合采用LSD方法在例8—1中，假定我们在分析之前就计划好要对品种1和品种3进行比较，看看这两个品种的产量之间是否有显著差异，这种情况下就适合采用LSD方法进行比较2018-9-25多重比较的LSD方法

2018-9-25多重比较的LSD方法

2018-9-25Fisher的LSD方法

(例题分析)library(agricolae)model_1w<-aov(产量~品种,data=example8_2)LSD<-LSD.test(model_1w,"品种");LSD

library(DescTools)PostHocTest(model_1w,method="lsd")2018-9-25Tukey-Kramer的HSD方法HSD是真实显著差异(honestlysignificantdifference)的缩写，因此也被称为真显著差异方法该检验方法是由JoneW.Tukey于1953年提出的，因此也被称为Tukey的HSD方法。由于Tukey的HSD方法要求各处理的样本量相同，当各处理的样本量不相同时，该方法就不再适用。20世纪50年代中期，C.Y.Kramer对Tukey的HSD方法做了一些修正，从而使其适用于样本量不同的情形。修正后的HSD检验称为Tukey-Kramer方法，简称为Tukey-Kramer的HSD方法该方法的适用场合是：研究者事先并未计划进行多重比较，只是在方差分析决绝原假设后，才需要对任意两个处理的均值进行比较，这时采用HSD方法比较合适2018-9-25Tukey-Kramer的HSD方法

2018-9-25Tukey-Kramer的HSD方法

2018-9-25Tukey-Kramer的HSD方法

(例题分析)#多重比较的TukeyHSD方法

#

多重比较的HSD方法（使用agricolae包输出其他信息）

#绘制配对差值置信区间的比较图

TukeyHSD(model_1w)library(agricolae)HSD<-HSD.test(model_1w,"品种");HSDplot(TukeyHSD(model_1w))8.3双因子方差分析

8.3.1数学模型

8.3.2主效应分析

8.3.3交互效应分析第8章方差分析2018-9-25双因子方差分析

(two-wayanalysisofvariance)

分析两个因子(因子A和因子B)对实验结果的影响如果两个因子对实验结果的影响是相互独立的，分别判断因子A和因子B对实验数据的单独影响，这时的双因子方差分析称为只考虑主效应的双因子方差分析或无重复双因子方差分析(Two-factorwithoutreplication)如果除了因子A和因子B对实验数据的单独影响外，两个因子的搭配还会对结果产生一种新的影响，这时的双因子方差分析称为考虑交互效应的双因子方差分析或可重复双因子方差分析

(Two-factorwithreplication)8.3.1数学模型8.3双因子方差分析2018-9-25双因子方差分析

(数学模型)设因子A有I种处理因子B有J种处理双因子方差分析可用下面的线性模型来表示

ij=0

8.3.2主效应分析

(maineffects)8.3双因子方差分析2018-9-25主效应分析

(效应检验)提出假设

2018-9-25主效应分析

(误差分解)2018-9-25主效应分析

(方差分析表)误差来源平方和SS自由度df均方MS检验统计量F因子A的处理效应SSA因子B的处理效应SSB误差SSE

总效应SST

2018-9-25主效应分析

(例题分析)【例8—5】(数据：example8_5.Rdata)假定在例8—1中，除了考虑品种对产量的影响外，还考虑施肥方式对产量的影响。假定有甲、乙两种施肥方式，这样3个小麦品种和两种施肥方式的搭配共有3×2=6种组合。如果选择30个地块进行实验，每一种搭配可以做5次实验，也就是每个品种（处理）的样本量为5，即相当于每个品种（处理）重复做了5次实验。实验取得的数据如表8—4所示。检验小麦品种和施肥方式对产量的影响是否显著(=0.05）2018-9-25将表8-4的短格式数据转为长格式数据

并另存为example8_5#加载数据table8_4<-read.csv("c:/example/ch8/table8_4.csv")table8_4<-cbind(table8_4,id=c(factor(1:10)))table8_4library(reshape)example8_5<-melt(table8_4,id.vars=c("id","施肥方式"))example8_5<-rename(example8_5,c(variable="品种",value="产量"))save(example8_5,file="C:/example/ch8/example8_5.RData")load("C:/example/ch8/example8_5.RData")example8_52018-9-25不同施肥方式下各品种产量的箱线图

(boxplot)load("C:/example/ch8/example8_5.RData");example8_5attach(example8_5)boxplot(产量~品种+施肥方式,col=c("gold","green","red"),ylab="产量",xlab="品种与施肥方式",data=example8_5)library(reshape)library(agricolae)mystats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))dfm<-melt(example8_5,measure.vars="产量",id.vars=c("品种","施肥方式"))cast(dfm,品种+施肥方式+variable~.,mystats)#加载数据#绘制品种和施肥方式的箱线图#按品种和施肥方式交叉分类计算均值和标准差2018-9-25#主效应方差分析结果

#主效应方差分析模型的参数估计方差分析表model_2wm<-aov(产量~品种+施肥方式)summary(model_2wm)model_2wm$coefficients2018-9-25主效应分析

(效应量分析)

2018-9-25主效应分析

(效应量分析)

2018-9-25主效应分析

(效应量分析)

2018-9-25#计算例8—5的效应量主效应分析

(效应量分析)model_2wm<-aov(产量~品种+施肥方式)library(DescTools)EtaSq(model_2wm,anova=T)8.3.3交互效应分析8.3双因子方差分析2018-9-25交互效应分析

(效应检验)提出假设

2018-9-25交互效应分析

(误差分解)2018-9-25交互效应分析

(方差分析表)误差来源平方和SS自由度df均方MS检验统计量F因子A的处理效应SSA因子B的处理效应SSBA、B的交互效应SSAB误差SSE

总效应SST

2018-9-25【例8-6】检验小麦品种、施肥方式及其交互效应对产量的影响是否显著#交互效应方差分析表#交互效应方差分析模型的参数估计

#绘制品种和施肥方式的主效应和交互效应图交互效应分析

(例题分析)attach(example8_5)fit<-aov(产量~品种+施肥方式+品种:施肥方式)summary(fit)fit$coefficientslibrary(HH)interaction2wt(产量~施肥方式+品种,data=example8_5)2018-9-25交互效应分析

(效应量分析)

2018-9-25#计算例8—6的效应量主效应分析

(效应量分析)model_2wi<-aov(产量~品种+施肥方式+品种:施肥方式,data=example8_5)library(DescTools)EtaSq(model_2wi,anova=T)2018-9-25主效应模型和交互效应模型的比较

(anova方法：例题8—4)#用anova函数比较模型model_2wm和model_2wimodel_2wm<-lm(产量~品种+施肥方式,data=example8_5)model_2wi<-lm(产量~品种+施肥方式+品种:施肥方式,data=example8_5)anova(model_2wm,model_2wi)8.4方差分析的假定及其检验

8.4.1正态性检验

8.4.2方差齐性检验第8章方差分析2018-9-25方差分析的基本假定正态性(normality)。每个总体都应服从正态分布，即对于因子的每一个水平，其观测值是来自正态分布总体的简单随机样本在例8—1中，要求每个品种的产量必须服从正态分布检验总体是否服从正态分布的方法有很多，包括对样本数据作直方图、茎叶图、箱线图、正态概率图做描述性判断，也可以进行非参数检验等方差齐性(homogeneityvariance)。各个总体的方差必须相同，对于分类变量的个水平，有

12=22=…=k2在例8—1中，要求不同品种的产量的方差都相同独立性(independence)。每个样本数据是来自因子各水平的独立样本(该假定不满足对结果影响较大)在例8—1中，3个样本数据是来自不同品种的3个独立样本8.4.1正态性检验8.4方差分析的假定及其检验2018-9-25正态性检验

(图示法)绘制因变量的正态概率图当每个处理的样本量足够大时，可以对每个样本绘制正态概率图来检查每个处理对应的总体是否服从正态分布当每个处理的样本量比较小时，正态概率图中的点很少，提供的正态性信息很有限。这时，可以将每个处理的样本数据合并绘制一个正态概率图来检验正态性2018-9-25正态性检验

(图示法)

#绘制每个品种产量数据的正态Q-Q图（数据：example8_1）load("C:/example/ch8/example8_1.RData")par(mfrow=c(1,3))qqnorm(example8_1$品种1,xlab="期望正态值",ylab="观察值",datax=TRUE,main="品种1的Q-Q图")qqline(example8_1$品种1,datax=TRUE)qqnorm(example8_1$品种2,xlab="期望正态值",ylab="观察值",datax=TRUE,main="品种2的Q-Q图")qqline(example8_1$品种2,datax=TRUE)qqnorm(example8_1$品种3,xlab="期望正态值",ylab="观察值",datax=TRUE,main="品种3的Q-Q图")qqline(example8_1$品种3,datax=TRUE)2018-9-25正态性检验

(图示法)

#绘制3个品种产量数据合并后的正态Q-Q图（数据：example8_2）load("C:/example/ch8/example8_2.RData")par(cex=.8,mai=c(.7,.7,.1,.1))qqnorm(example8_2$产量,xlab="期望正态值",ylab="观察值",data=TRUE,main="")qqline(example8_2$产量,datax=TRUE,col="red",lwd=2)op<-par(fig=c(.08,.5,.5,.98),new=TRUE)hist(example8_2$产量,xlab="产量",ylab="",freq=FALSE,col="lightblue",cex.axis=0.7,cex.lab=0.7,main="")lines(density(example8_2$产量),col="red",lwd=2)box()2018-9-25正态性检验

(检验法)当样本量较小时，正态概率图的应用就会受到很大限制，这时可以使用标准的统计检验如Shapiro—Wilk检验、Kolmogorov-Smirnov检验等，均可以做正态性检验。这些检验的原假设是因变量服从正态分布如果检验获得的P值小于指定的显著性水平，则拒绝原假设，表明总体不服从正态分布，如果P值较大不能拒绝原假设时，可以认为总体满足正态分布这些检验对正态性的轻微偏离是敏感的，检验往往导致拒绝原假设。而方差分析对正态性的要求则相对比较宽松，当正态性略微不满足时，对分析结果的影响不是很大。因此，实际中应谨慎使用这些检验2018-9-25正态性检验

(检验法)

##每个品种产量的正态性检验#品种1的正态性检验：#品种2的正态性检验：#品种3的正态性检验：load("C:/example/ch8/example8_2.RData")attach(example8_2)shapiro.test(产量[品种=='品种1'])ks.test(产量[品种=='品种1'],"pnorm",mean(产量[品种=='品种1']),sd(产量[品种=='品种1']))shapiro.test(产量[品种=='品种2'])ks.test(产量[品种=='品种2'],"pnorm",mean(产量[品种=='品种2']),sd(产量[品种=='品种2']))shapiro.test(产量[品种=='品种3'])ks.test(产量[品种=='品种3'],"pnorm",mean(产量[品种=='品种3']),sd(产量[品种=='品种3']))2018-9-25正态性检验

(检验法)

##三个品种产量数据合并后的检验#shapiro正态性检验#K-S正态性检验

load("C:/example/ch8/example8_2.RData")attach(example8_2)shapiro.test(产量)ks.test(产量,"pnorm",mean(产量),sd(产量))2018-9-25正态性检验

(检验法)

##

双因子方差分析：Shapiro-Wilk正态性检验和k-s检验（数据：example8_5）#施肥方式甲的产量正态性检验：#施肥方式乙的产量正态性检验：load("C:/example/ch8/example8_5.RData")attach(example8_5)shapiro.test(产量[施肥方式=='甲'])ks.test(产量[施肥方式=='甲'],"pnorm",mean(产量[施肥方式=='甲']),sd(产量[施肥方式=='甲']))shapiro.test(产量[施肥方式=='乙'])ks.test(产量[施肥方式=='乙'],"pnorm",mean(产量[施肥方式=='乙']),sd(产量[施肥方式=='乙']))8.4.2方差齐性检验8.4方差分析的假定及其检验2018-9-25方差齐性检验方差齐性(homogeneityvariance)。假定各个总体的方差必须相同，即：

12=22=…=I2在例8—1中，要求不同品种的产量的方差都相同检验方差齐性方法：图示法和检验法

检验方差齐性的图形有箱线图和残差图等Bartlett方差齐性检验和Levene方差齐性检验

2018-9-25方差齐性检验

(图示法)#绘制例8—2方差分析的残差图和残差的Q-Q（资料：example8_2）

#绘制例8—6方差分析的残差图和残差的Q-Q（资料：example8_5）load("C:/example/ch8/example8_2.RData")model_1w<-aov(产量~品种,data=example8_2)par(mfrow=c(1,2),mai=c(0.5,0.5,0.2,0.1),cex=0.6,cex.main=0.7)

plot(model_1w,which=1:2)load("C:/example/ch8/example8_5.RData")model_2wi<-aov(产量~品种+施肥方式+品种:施肥方式,data=example8_5)par(m