《統(tǒng)計學(xué)-基于R》(08)第8章-方差分析(R3)

數(shù)據(jù)分析

(方法與案例)

作者賈俊平版權(quán)所有違者必究StatisticswithR統(tǒng)計學(xué)R語言第8章方差分析8.1

方差分析的原理8.2單因子方差分析8.3雙因子方差分析8.4方差分析的假定及其檢驗8.5單因子方差分析的非參數(shù)方法NOVAA8.1方差分析的原理

8.1.1什么是方差分析

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

(analysisofvariance)

方差分析的基本原理是在20世紀(jì)20年代由英國統(tǒng)計學(xué)家RonaldA.Fisher在進行實驗設(shè)計時為解釋實驗數(shù)據(jù)而首先引入的分析各分類自變量對數(shù)值因變量影響的一種統(tǒng)計方法研究分類型自變量對數(shù)值型因變量的影響

一個或多個分類型自變量兩個或多個(k個)處理水平或分類一個數(shù)值型因變量有單因子方差分析和雙因子方差分析單因子方差分析：涉及一個分類的自變量雙因子方差分析：涉及兩個分類的自變量2018-9-25什么是方差分析

(例題分析)品種1品種2品種3817176827279797277816676787278897789928187877784857387867987【例8-1】（數(shù)據(jù)：example8_1.RData）為分析小麥品種對產(chǎn)量的影響，一家研究機構(gòu)挑選了3個小麥品種：品種1、品種2、品種3，然后選擇條件和面積相同的30個地塊，每個品種在10個地塊上試種，實驗獲得的產(chǎn)量數(shù)據(jù)如表8-1所示。2018-9-25什么是方差分析

(例題分析)分析“小麥品種”對“產(chǎn)量”的影響如果只分析品種一個因子對產(chǎn)量的影響，則稱為單因子方差分析(one-wayanalysisofvariance)如果兩個因子對產(chǎn)量的單獨影響，但不考慮它們對產(chǎn)量的交互效應(yīng)(interaction)，則稱為只考慮主效應(yīng)(maineffect)的雙因子方差分析，或稱為無重復(fù)雙因子分析(two-factorwithoutreplication)如果除了考慮兩個因子對產(chǎn)量的單獨影響外，還考慮二者對產(chǎn)量的交互效應(yīng)，則稱為考慮交互效應(yīng)的雙因子方差分析，或稱為可重復(fù)雙因子分析(two-factorwithreplication)8.1.2誤差分解8.1方差分析的原理2018-9-25方差分析的基本原理

(誤差分解)總誤差(totalerror)反映全部觀測數(shù)據(jù)的誤差所抽取的全部30個地塊的產(chǎn)量之間差異處理誤差(treatmenterror)—組間誤差(between-grouperror)由于不同處理造成的誤差，它反映了處理(品種)對觀測數(shù)據(jù)(產(chǎn)量)的影響，因此稱為處理效應(yīng)(treatmenteffect)隨機誤差(randomerror)—組內(nèi)誤差(within-grouperror)由于隨機因子造成的誤差，也簡稱為誤差(error)2018-9-25方差分析的基本原理

(誤差分解)數(shù)據(jù)的誤差用平方和(sumofsquares)表示，記為SS總平方和(sumofsquaresfortotal)，記為SST反映全部數(shù)據(jù)總誤差大小的平方和抽取的全部30個地塊產(chǎn)量之間的誤差平方和處理平方和(treatmentsumofsquares)，記為SSA反映處理誤差大小的平方和也稱為組間平方和(between-groupsumofsquares)誤差平方和(sumofsquaresoferror)，記為SSE反映隨機誤差大小的平方和稱為誤差平方和也稱為組內(nèi)平方和(within-groupsumofsquares)2018-9-25方差分析的基本原理

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

(誤差分析)方差分析的基本原理就是要分析數(shù)據(jù)的總誤差中有沒有處理誤差。如果處理(超市的不同位置)對觀測數(shù)據(jù)(銷售額)沒有顯著影響，意味著沒有處理誤差。這時，每種處理所對應(yīng)的總體均值(

i)應(yīng)該相等如果存在處理誤差，每種處理所對應(yīng)的總體均值(

i)至少有一對不相等就例8—1而言，在只考慮品種一個因子的情況下，方差分析也就是要檢驗下面的假設(shè)H0：

1

2

3

H1：

1,

2,

3

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

8.2.1數(shù)學(xué)模型

8.2.2效應(yīng)檢驗

8.2.3效應(yīng)量分析

8.2.4多重比較第8章方差分析8.2.1數(shù)學(xué)模型8.2單因子方差分析2018-9-25單因子方差分析

(數(shù)學(xué)模型)設(shè)因子A有I種處理(比如品種有“品種1”、“品種2”、“品種3”3種處理)，單因子方差分析可用下面的線性模型來表示

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

(數(shù)學(xué)模型)設(shè)全部觀測數(shù)據(jù)的總均值為

，第i個處理效應(yīng)用第i個處理均值與總均值的差(

i-)

表示，記為

i，即i=i-

。這樣，第i個處理均值被分解成i=i+，方差分析模型可以表達為

8.2.2效應(yīng)檢驗8.2單因子方差分析2018-9-25提出假設(shè)一般提法H0

：

i

=0

(i=1,2,…,I)沒有處理效應(yīng)

H1：

i

至少有一個不等于0有處理效應(yīng)

注意：拒絕原假設(shè)，只表明至少有一個處理的效應(yīng)顯著，并不意味著所有的粗粒的效應(yīng)都顯著2018-9-25構(gòu)造檢驗的統(tǒng)計量F單因子方差分析的方差分析表誤差來源平方和SS自由度df均方MS檢驗統(tǒng)計量F處理效應(yīng)誤差

總效應(yīng)

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

(例題分析)

#將表8-1的短格式數(shù)據(jù)轉(zhuǎn)為長格式數(shù)據(jù)，并顯示前6列l(wèi)oad("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="產(chǎn)量"))save(example8_2,file="C:/example/ch8/example8_2.RData")example8_22018-9-25#繪制3個品種數(shù)據(jù)產(chǎn)量的箱線圖

#計算描述統(tǒng)計量3個品種產(chǎn)量的描述性分析load("C:/example/ch8/example8_2.RData")attach(example8_2)boxplot(產(chǎn)量~品種,data=example8_2,col="gold",main="",ylab="產(chǎn)量",xlab="品種")my_summary<-function(x){with(x,data.frame("均值"=mean(產(chǎn)量),"標(biāo)準(zhǔn)差"=sd(產(chǎn)量),n=length(產(chǎn)量)))}library(plyr)ddply(example8_2,.(品種),my_summary)2018-9-25#方差分析表#方差分析模型的參數(shù)估計#繪制均值圖方差分析表和參數(shù)估計attach(example8_2)model_1w<-aov(產(chǎn)量~品種)summary(model_1w)model_1w$coefficientslibrary(gplots)plotmeans(產(chǎn)量~品種,data=example8_2)8.2.3效應(yīng)量分析8.2單因子方差分析2018-9-25效應(yīng)量分析

2018-9-25library(DescTools)

model_1w<-aov(產(chǎn)量~品種)

EtaSq(model_1w,anova=T)效應(yīng)量分析

【例8-2】的效應(yīng)量8.2.3多重比較8.2單因子方差分析2018-9-25多重比較的意義通過對總體均值之間的配對比較來進一步檢驗到底哪些均值之間存在差異比較方法有多種，如Fisher的LSD方法、Tukey-Kramer的HSD方法等2018-9-25Fisher的LSD方法LSD是最小顯著差異(leastsignificantdifference)的縮寫，該檢驗方法是由統(tǒng)計學(xué)家Fisher提出來的，因此也稱為Fisher的最小顯著差異方法，簡稱LSD方法LSD的適用場合：如果研究者在事先就已經(jīng)計劃好要對某對或某幾對均值進行比較，不管方差分析的結(jié)果如何(拒絕或不拒絕原假設(shè))，都要進行比較，這時適合采用LSD方法在例8—1中，假定我們在分析之前就計劃好要對品種1和品種3進行比較，看看這兩個品種的產(chǎn)量之間是否有顯著差異，這種情況下就適合采用LSD方法進行比較2018-9-25多重比較的LSD方法

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

2018-9-25Fisher的LSD方法

(例題分析)library(agricolae)model_1w<-aov(產(chǎn)量~品種,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方法要求各處理的樣本量相同，當(dāng)各處理的樣本量不相同時，該方法就不再適用。20世紀(jì)50年代中期，C.Y.Kramer對Tukey的HSD方法做了一些修正，從而使其適用于樣本量不同的情形。修正后的HSD檢驗稱為Tukey-Kramer方法，簡稱為Tukey-Kramer的HSD方法該方法的適用場合是：研究者事先并未計劃進行多重比較，只是在方差分析決絕原假設(shè)后，才需要對任意兩個處理的均值進行比較，這時采用HSD方法比較合適2018-9-25Tukey-Kramer的HSD方法

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

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

(例題分析)#多重比較的TukeyHSD方法

#

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

#繪制配對差值置信區(qū)間的比較圖

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

8.3.1數(shù)學(xué)模型

8.3.2主效應(yīng)分析

8.3.3交互效應(yīng)分析第8章方差分析2018-9-25雙因子方差分析

(two-wayanalysisofvariance)

分析兩個因子(因子A和因子B)對實驗結(jié)果的影響如果兩個因子對實驗結(jié)果的影響是相互獨立的，分別判斷因子A和因子B對實驗數(shù)據(jù)的單獨影響，這時的雙因子方差分析稱為只考慮主效應(yīng)的雙因子方差分析或無重復(fù)雙因子方差分析(Two-factorwithoutreplication)如果除了因子A和因子B對實驗數(shù)據(jù)的單獨影響外，兩個因子的搭配還會對結(jié)果產(chǎn)生一種新的影響，這時的雙因子方差分析稱為考慮交互效應(yīng)的雙因子方差分析或可重復(fù)雙因子方差分析

(Two-factorwithreplication)8.3.1數(shù)學(xué)模型8.3雙因子方差分析2018-9-25雙因子方差分析

(數(shù)學(xué)模型)設(shè)因子A有I種處理因子B有J種處理雙因子方差分析可用下面的線性模型來表示

ij=0

8.3.2主效應(yīng)分析

(maineffects)8.3雙因子方差分析2018-9-25主效應(yīng)分析

(效應(yīng)檢驗)提出假設(shè)

2018-9-25主效應(yīng)分析

(誤差分解)2018-9-25主效應(yīng)分析

(方差分析表)誤差來源平方和SS自由度df均方MS檢驗統(tǒng)計量F因子A的處理效應(yīng)SSA因子B的處理效應(yīng)SSB誤差SSE

總效應(yīng)SST

2018-9-25主效應(yīng)分析

(例題分析)【例8—5】(數(shù)據(jù)：example8_5.Rdata)假定在例8—1中，除了考慮品種對產(chǎn)量的影響外，還考慮施肥方式對產(chǎn)量的影響。假定有甲、乙兩種施肥方式，這樣3個小麥品種和兩種施肥方式的搭配共有3×2=6種組合。如果選擇30個地塊進行實驗，每一種搭配可以做5次實驗，也就是每個品種（處理）的樣本量為5，即相當(dāng)于每個品種（處理）重復(fù)做了5次實驗。實驗取得的數(shù)據(jù)如表8—4所示。檢驗小麥品種和施肥方式對產(chǎn)量的影響是否顯著(=0.05）2018-9-25將表8-4的短格式數(shù)據(jù)轉(zhuǎn)為長格式數(shù)據(jù)

并另存為example8_5#加載數(shù)據(jù)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="產(chǎn)量"))save(example8_5,file="C:/example/ch8/example8_5.RData")load("C:/example/ch8/example8_5.RData")example8_52018-9-25不同施肥方式下各品種產(chǎn)量的箱線圖

(boxplot)load("C:/example/ch8/example8_5.RData");example8_5attach(example8_5)boxplot(產(chǎn)量~品種+施肥方式,col=c("gold","green","red"),ylab="產(chǎn)量",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="產(chǎn)量",id.vars=c("品種","施肥方式"))cast(dfm,品種+施肥方式+variable~.,mystats)#加載數(shù)據(jù)#繪制品種和施肥方式的箱線圖#按品種和施肥方式交叉分類計算均值和標(biāo)準(zhǔn)差2018-9-25#主效應(yīng)方差分析結(jié)果

#主效應(yīng)方差分析模型的參數(shù)估計方差分析表model_2wm<-aov(產(chǎn)量~品種+施肥方式)summary(model_2wm)model_2wm$coefficients2018-9-25主效應(yīng)分析

(效應(yīng)量分析)

2018-9-25主效應(yīng)分析

(效應(yīng)量分析)

2018-9-25主效應(yīng)分析

(效應(yīng)量分析)

2018-9-25#計算例8—5的效應(yīng)量主效應(yīng)分析

(效應(yīng)量分析)model_2wm<-aov(產(chǎn)量~品種+施肥方式)library(DescTools)EtaSq(model_2wm,anova=T)8.3.3交互效應(yīng)分析8.3雙因子方差分析2018-9-25交互效應(yīng)分析

(效應(yīng)檢驗)提出假設(shè)

2018-9-25交互效應(yīng)分析

(誤差分解)2018-9-25交互效應(yīng)分析

(方差分析表)誤差來源平方和SS自由度df均方MS檢驗統(tǒng)計量F因子A的處理效應(yīng)SSA因子B的處理效應(yīng)SSBA、B的交互效應(yīng)SSAB誤差SSE

總效應(yīng)SST

2018-9-25【例8-6】檢驗小麥品種、施肥方式及其交互效應(yīng)對產(chǎn)量的影響是否顯著#交互效應(yīng)方差分析表#交互效應(yīng)方差分析模型的參數(shù)估計

#繪制品種和施肥方式的主效應(yīng)和交互效應(yīng)圖交互效應(yīng)分析

(例題分析)attach(example8_5)fit<-aov(產(chǎn)量~品種+施肥方式+品種:施肥方式)summary(fit)fit$coefficientslibrary(HH)interaction2wt(產(chǎn)量~施肥方式+品種,data=example8_5)2018-9-25交互效應(yīng)分析

(效應(yīng)量分析)

2018-9-25#計算例8—6的效應(yīng)量主效應(yīng)分析

(效應(yīng)量分析)model_2wi<-aov(產(chǎn)量~品種+施肥方式+品種:施肥方式,data=example8_5)library(DescTools)EtaSq(model_2wi,anova=T)2018-9-25主效應(yīng)模型和交互效應(yīng)模型的比較

(anova方法：例題8—4)#用anova函數(shù)比較模型model_2wm和model_2wimodel_2wm<-lm(產(chǎn)量~品種+施肥方式,data=example8_5)model_2wi<-lm(產(chǎn)量~品種+施肥方式+品種:施肥方式,data=example8_5)anova(model_2wm,model_2wi)8.4方差分析的假定及其檢驗

8.4.1正態(tài)性檢驗

8.4.2方差齊性檢驗第8章方差分析2018-9-25方差分析的基本假定正態(tài)性(normality)。每個總體都應(yīng)服從正態(tài)分布，即對于因子的每一個水平，其觀測值是來自正態(tài)分布總體的簡單隨機樣本在例8—1中，要求每個品種的產(chǎn)量必須服從正態(tài)分布檢驗總體是否服從正態(tài)分布的方法有很多，包括對樣本數(shù)據(jù)作直方圖、莖葉圖、箱線圖、正態(tài)概率圖做描述性判斷，也可以進行非參數(shù)檢驗等方差齊性(homogeneityvariance)。各個總體的方差必須相同，對于分類變量的個水平，有

12=22=…=k2在例8—1中，要求不同品種的產(chǎn)量的方差都相同獨立性(independence)。每個樣本數(shù)據(jù)是來自因子各水平的獨立樣本(該假定不滿足對結(jié)果影響較大)在例8—1中，3個樣本數(shù)據(jù)是來自不同品種的3個獨立樣本8.4.1正態(tài)性檢驗8.4方差分析的假定及其檢驗2018-9-25正態(tài)性檢驗

(圖示法)繪制因變量的正態(tài)概率圖當(dāng)每個處理的樣本量足夠大時，可以對每個樣本繪制正態(tài)概率圖來檢查每個處理對應(yīng)的總體是否服從正態(tài)分布當(dāng)每個處理的樣本量比較小時，正態(tài)概率圖中的點很少，提供的正態(tài)性信息很有限。這時，可以將每個處理的樣本數(shù)據(jù)合并繪制一個正態(tài)概率圖來檢驗正態(tài)性2018-9-25正態(tài)性檢驗

(圖示法)

#繪制每個品種產(chǎn)量數(shù)據(jù)的正態(tài)Q-Q圖（數(shù)據(jù)：example8_1）load("C:/example/ch8/example8_1.RData")par(mfrow=c(1,3))qqnorm(example8_1$品種1,xlab="期望正態(tài)值",ylab="觀察值",datax=TRUE,main="品種1的Q-Q圖")qqline(example8_1$品種1,datax=TRUE)qqnorm(example8_1$品種2,xlab="期望正態(tài)值",ylab="觀察值",datax=TRUE,main="品種2的Q-Q圖")qqline(example8_1$品種2,datax=TRUE)qqnorm(example8_1$品種3,xlab="期望正態(tài)值",ylab="觀察值",datax=TRUE,main="品種3的Q-Q圖")qqline(example8_1$品種3,datax=TRUE)2018-9-25正態(tài)性檢驗

(圖示法)

#繪制3個品種產(chǎn)量數(shù)據(jù)合并后的正態(tài)Q-Q圖（數(shù)據(jù)：example8_2）load("C:/example/ch8/example8_2.RData")par(cex=.8,mai=c(.7,.7,.1,.1))qqnorm(example8_2$產(chǎn)量,xlab="期望正態(tài)值",ylab="觀察值",data=TRUE,main="")qqline(example8_2$產(chǎn)量,datax=TRUE,col="red",lwd=2)op<-par(fig=c(.08,.5,.5,.98),new=TRUE)hist(example8_2$產(chǎn)量,xlab="產(chǎn)量",ylab="",freq=FALSE,col="lightblue",cex.axis=0.7,cex.lab=0.7,main="")lines(density(example8_2$產(chǎn)量),col="red",lwd=2)box()2018-9-25正態(tài)性檢驗

(檢驗法)當(dāng)樣本量較小時，正態(tài)概率圖的應(yīng)用就會受到很大限制，這時可以使用標(biāo)準(zhǔn)的統(tǒng)計檢驗如Shapiro—Wilk檢驗、Kolmogorov-Smirnov檢驗等，均可以做正態(tài)性檢驗。這些檢驗的原假設(shè)是因變量服從正態(tài)分布如果檢驗獲得的P值小于指定的顯著性水平，則拒絕原假設(shè)，表明總體不服從正態(tài)分布，如果P值較大不能拒絕原假設(shè)時，可以認(rèn)為總體滿足正態(tài)分布這些檢驗對正態(tài)性的輕微偏離是敏感的，檢驗往往導(dǎo)致拒絕原假設(shè)。而方差分析對正態(tài)性的要求則相對比較寬松，當(dāng)正態(tài)性略微不滿足時，對分析結(jié)果的影響不是很大。因此，實際中應(yīng)謹(jǐn)慎使用這些檢驗2018-9-25正態(tài)性檢驗

(檢驗法)

##每個品種產(chǎn)量的正態(tài)性檢驗#品種1的正態(tài)性檢驗：#品種2的正態(tài)性檢驗：#品種3的正態(tài)性檢驗：load("C:/example/ch8/example8_2.RData")attach(example8_2)shapiro.test(產(chǎn)量[品種=='品種1'])ks.test(產(chǎn)量[品種=='品種1'],"pnorm",mean(產(chǎn)量[品種=='品種1']),sd(產(chǎn)量[品種=='品種1']))shapiro.test(產(chǎn)量[品種=='品種2'])ks.test(產(chǎn)量[品種=='品種2'],"pnorm",mean(產(chǎn)量[品種=='品種2']),sd(產(chǎn)量[品種=='品種2']))shapiro.test(產(chǎn)量[品種=='品種3'])ks.test(產(chǎn)量[品種=='品種3'],"pnorm",mean(產(chǎn)量[品種=='品種3']),sd(產(chǎn)量[品種=='品種3']))2018-9-25正態(tài)性檢驗

(檢驗法)

##三個品種產(chǎn)量數(shù)據(jù)合并后的檢驗#shapiro正態(tài)性檢驗#K-S正態(tài)性檢驗

load("C:/example/ch8/example8_2.RData")attach(example8_2)shapiro.test(產(chǎn)量)ks.test(產(chǎn)量,"pnorm",mean(產(chǎn)量),sd(產(chǎn)量))2018-9-25正態(tài)性檢驗

(檢驗法)

##

雙因子方差分析：Shapiro-Wilk正態(tài)性檢驗和k-s檢驗（數(shù)據(jù)：example8_5）#施肥方式甲的產(chǎn)量正態(tài)性檢驗：#施肥方式乙的產(chǎn)量正態(tài)性檢驗：load("C:/example/ch8/example8_5.RData")attach(example8_5)shapiro.test(產(chǎn)量[施肥方式=='甲'])ks.test(產(chǎn)量[施肥方式=='甲'],"pnorm",mean(產(chǎn)量[施肥方式=='甲']),sd(產(chǎn)量[施肥方式=='甲']))shapiro.test(產(chǎn)量[施肥方式=='乙'])ks.test(產(chǎn)量[施肥方式=='乙'],"pnorm",mean(產(chǎn)量[施肥方式=='乙']),sd(產(chǎn)量[施肥方式=='乙']))8.4.2方差齊性檢驗8.4方差分析的假定及其檢驗2018-9-25方差齊性檢驗方差齊性(homogeneityvariance)。假定各個總體的方差必須相同，即：

12=22=…=I2在例8—1中，要求不同品種的產(chǎn)量的方差都相同檢驗方差齊性方法：圖示法和檢驗法

檢驗方差齊性的圖形有箱線圖和殘差圖等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(產(chǎn)量~品種,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(產(chǎn)量~品種+施肥方式+品種:施肥方式,data=example8_5)par(m