《統(tǒng)計(jì)計(jì)算及統(tǒng)計(jì)軟件》課件ch6

問題1：某個農(nóng)科所要研究某農(nóng)作物產(chǎn)量，考慮了四個因素：品種、密度、肥料、土壤；為方便分析，每個因素設(shè)定了三個不同的類別(或水平)試回答：(1)全面完成這項(xiàng)工作需要安排幾次試驗(yàn)？(2)如何安排試驗(yàn)?zāi)軌蚬?jié)省成本又不影響效果？方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）問題2：某個教師的授課班中有四個專業(yè)，某次測試的抽樣成績?nèi)缬冶恚瑔枺?1)四個專業(yè)的學(xué)生學(xué)習(xí)效果是否無異？(2)若老師有多個，那老師、專業(yè)分別有差異否？方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）第六章方差分析與正交試驗(yàn)設(shè)計(jì)

概述

方差分析6.1單因素方差分析6.2雙因素方差分析

分析前提6.3方差齊性與均值差異的檢驗(yàn)

正交試驗(yàn)6.4正交試驗(yàn)設(shè)計(jì)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）概述

方差分析與正交試驗(yàn)源自農(nóng)業(yè)生產(chǎn)，由R.A.Fisher于上世紀(jì)20年代創(chuàng)立，目前已得到充分的推廣和應(yīng)用。

方差分析的前提基于正態(tài)分布和同方差假定，研究的是因素(自變量)與試驗(yàn)指標(biāo)(因變量)之間的關(guān)系，確定哪些因素對試驗(yàn)指標(biāo)有顯著的影響。

正交試驗(yàn)通過巧妙地安排試驗(yàn)方案和次數(shù)，在降低成本的同時確定較優(yōu)(最優(yōu))的試驗(yàn)方案。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）概述

依據(jù)影響試驗(yàn)指標(biāo)的因素的個數(shù)，將方差分析區(qū)分為單因素、雙因素和多因素。通常用大寫字母表示因素，如A,B,C,…，而用X表示實(shí)驗(yàn)指標(biāo)。

每個因素根據(jù)需要劃分成不同的狀態(tài)或等級或?qū)哟?，稱為水平，分別記為A1,A2,…,Ar。

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）概述

試驗(yàn)結(jié)果的不同可能是由不同因素的不同水平造成的(因素效應(yīng))，也可能是由系統(tǒng)性誤差(試驗(yàn)誤差)，方差分析的主要任務(wù)就是分離因素效應(yīng)和試驗(yàn)誤差，確定各因素的主次順序和影響程度。

實(shí)際操作中，不同水平下所作試驗(yàn)的次數(shù)不盡相同，若相同則稱為等重復(fù)試驗(yàn)，否則為不等重復(fù)試驗(yàn)。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析例6.1.1肥料種類優(yōu)選試驗(yàn)。某農(nóng)科所為比較三種不同肥料對水稻產(chǎn)量的影響，進(jìn)行如下試驗(yàn)：選一塊肥沃程度較均勻的土地，將其分割成9塊，為減少肥沃程度影響，按照表6.1.3所示安排試驗(yàn)。推斷肥料種類對水稻產(chǎn)量具有顯著的影響？(取顯著性水平)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析基本假設(shè)： 1.每個水平下的試驗(yàn)指標(biāo)看作一個總體，且服從正態(tài)分布，試驗(yàn)結(jié)果看成總體的一個抽樣，且各水平試驗(yàn)相互獨(dú)立。 2.每個水平對應(yīng)總體的方差相等(方差齊性)

則上述問題變成檢驗(yàn)各總體的均值是否相等，即方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析

顯然，若拒絕Ho，則說明不同水平的效應(yīng)之間存在顯著差異。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析

如何從數(shù)據(jù)中分離出因素效應(yīng)和誤差效應(yīng)是執(zhí)行檢驗(yàn)的關(guān)鍵，該技術(shù)依賴于正態(tài)前提下的總離差平方和分解及自由度分解理論。

平方和分解：

交互項(xiàng)為0SA和SE相互獨(dú)立，且方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）參考：孫榮恒.應(yīng)用數(shù)理統(tǒng)計(jì)[M].科學(xué)出版社,2003,P18-216.1單因素方差分析SA和SE相互獨(dú)立，且則可構(gòu)造檢驗(yàn)統(tǒng)計(jì)量

若Ho成立，則各水平無差異，則組間的偏差很小，偏差主要為組內(nèi)偏差，即F值應(yīng)很小。

顯然若F值越大說明因素的效應(yīng)越明顯，即偏差是由因素的水平不同造成的，從而拒絕原假設(shè)Ho。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析#H0成立時的卡方分布模擬，N(10,1)為原分布，水平=5，試驗(yàn)次數(shù)均等為100getchisq=function(n=500,ni=100,r=5,mu=10,xigma=1){SE=numeric(n);SA=numeric(n);mat=matrix(0,ni,r)for(iin1:n){for(jin1:r)mat[,j]=rnorm(ni,mu,xigma)mx=apply(mat,2,mean);ax=mean(mat);md=apply(mat,2,var)*(ni-1)SE[i]=sum(md);SA[i]=sum((mx-ax)^2*(ni-1))}ST=SE+SAp1=ks.test(ST,'pchisq',n-1)$p.value;p2=ks.test(SA,'pchisq',r-1)$p.value;p3=ks.test(SE,'pchisq',n-r)$p.valuereturn(c(p1,p2,p3))}allok=function(x,level=0.1)return(sum(x>level))a=rep(500,100);b=sapply(a,getchisq)sum(apply(t(b),1,allok)==3)/100方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析在R中可通過aov函數(shù)實(shí)現(xiàn)方差分析.針對例題數(shù)據(jù)，R程序如下：aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)dat1=data.frame(x=c(94,91,75,62,68,50,78,65,80),A=factor(rep(1:3,each=3)))aov.out1=aov(x~A,data=dat1)#可以查閱names(aov.dat)summary(aov.out1)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.1單因素方差分析aov1.table=function(aov.out,alpha=0.05){a=summary(aov.out)#得到方差分析結(jié)果b=a[[1]];cols=ncol(b);#去掉列表中的數(shù)據(jù)框addrow=c(sum(b[,1]),sum(b[,2]),sum(b[,2])/sum(b[,1]),qf(1-alpha,b[1,1],b[2,1]),alpha)a[[1]]['Total',]=addrow#增補(bǔ)行(自由度,總和,均方和,F分位點(diǎn),顯著性水平)return(a[[1]][,c(2,1,3:cols)])#返回方差分析結(jié)果,前兩列交換了}#下面一句通過因子確定數(shù)據(jù)歸屬aov.out1=aov(x~A,data=dat1)#得到方差分析對象aov1.table(aov.out1)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2雙因素方差分析

因素A有r個水平,因素B有s個水平，每個交互項(xiàng)AiBj看作一個總體，分別作m次試驗(yàn)，得到結(jié)果如下：方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2雙因素方差分析

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）AiBj總效應(yīng)減去Ai效應(yīng)和Bj效應(yīng)得到AiBj的交互效應(yīng)。6.2雙因素方差分析

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2雙因素方差分析

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）課堂練習(xí)：

數(shù)據(jù)1：

數(shù)據(jù)2：方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）P(F(3,6)>9.7795)=0.01,P(F(2,6)>10.9247)=0.01)你的結(jié)論：P(F(2,12)>6.926)=0.01,P(F(6,12)>2.33)=0.1)你的結(jié)論：dat1=c(12,13,13,14,12,15,14,15,13,15,16,16)A=gl(3,4,12)B=gl(4,1,12)summary(aov(dat1~A+B))dat2=c(12,13,13,14,13,14,14,14,12,14,15,14,14,16,15,16,13,15,15,16,16,17,16,15)A=gl(3,8,24)B=gl(4,2,24)summary(aov(dat2~A+B+A:B))方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.3方差齊性檢驗(yàn)

方差齊性檢驗(yàn)就是檢驗(yàn)數(shù)據(jù)在不同因素(不同水平)下的方差是否相同，其中最常用的方法是Bartlett檢驗(yàn)。該方法對分布的非正態(tài)性比較敏感，只要總體方差有顯著差別或者總體分布的偏斜程度有所不同，Bartlett檢驗(yàn)的結(jié)果都可能顯著。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.3方差齊性檢驗(yàn)

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）在R中，可以使用函數(shù)bartlett.test實(shí)現(xiàn)方差齊性檢驗(yàn)，其定義如下：bartlett.test(x,g,subset,...)

bartlett.test(formula,data,...)

顯然統(tǒng)計(jì)量值越大越應(yīng)該拒絕6.3方差齊性檢驗(yàn)經(jīng)以下代碼模擬，發(fā)現(xiàn)這個近似分布還比較靠譜。n=300;m=6;k=100;c=1+(m/(n/m-1)-1/(n-m))/(3*(m-1))chi=numeric(k)for(iin1:k){mat=matrix(rnorm(n),nrow=n/m);s2=apply(mat,2,var)ni=rep(n/m,m)c2=sum((ni-1)*log(s2));c1=(n-m)*log(sum((ni-1)*s2)/(n-m))chi[i]=(c1-c2)/c}ks.test(chi,'pchisq',m-1)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.4正交試驗(yàn)設(shè)計(jì)（orthogonaldesign） 簡稱正交設(shè)計(jì)（orthoplan），它研究的是如何利用正交表(orthogonaltable)工具來安排試驗(yàn)，用盡可能少的試驗(yàn)次數(shù)獲得盡可能多的典型數(shù)據(jù)，然后進(jìn)行特定的統(tǒng)計(jì)分析，最終得出最優(yōu)或較優(yōu)的試驗(yàn)條件或方案主要內(nèi)容：

1.正交表

2.試驗(yàn)安排

3.直觀分析、極差分析和方差分析方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.4.1正交表

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）正交表是一種特殊的二維表格,通常表示成

n表示試驗(yàn)次數(shù),m表示因素(因子)個數(shù),r表示水平數(shù)最多安排m個因子，每個因子為r個水平的試驗(yàn)2.表示最多安排n行m列，表中數(shù)據(jù)都是1,…,r的數(shù)字3.全部試驗(yàn)是r^m個，但是試驗(yàn)只安排n個全面試驗(yàn)81個，安排9個試驗(yàn)，最多支持4因子3水平。6.4.1正交表

方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）觀察表格中的數(shù)字(特點(diǎn))：

1.每列數(shù)字出現(xiàn)的次數(shù)相等

2.行可交換，列可交換

3.任意兩列構(gòu)成的數(shù)對出現(xiàn)的次數(shù)相等，且搭配均衡全面。6.4.2正交試驗(yàn)基本步驟方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）

1.由因子個數(shù)和水平，確定正交表，原則是水平數(shù)因子數(shù)(找出滿足要求的因子數(shù)最小的)正交表

r=ro,m>=m0,min{m}2.表頭設(shè)計(jì)，將各因素放置在不同的列(是否隨機(jī)）可能出現(xiàn)空列(此時一般將其作為誤差列)。

3.嚴(yán)格按照實(shí)驗(yàn)條件(因素水平組合)進(jìn)行試驗(yàn)，記錄試驗(yàn)結(jié)果。

4.進(jìn)行直觀分析、極差分析和方差分析確定最佳試驗(yàn)條件和因素的影響次序6.4.2正交試驗(yàn)例子分析每個因素有三個水平應(yīng)選取正交表________(溫度，時間，紙漿溶度)方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.4.2正交試驗(yàn)例子分析—直觀分析方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）所謂直接分析：就是從正交試驗(yàn)中找出試驗(yàn)結(jié)果最好的那一組試驗(yàn)及其條件。特別注意的是：試驗(yàn)結(jié)果規(guī)定為越大越好（以大為好)。6.4.2正交試驗(yàn)例子分析—極差分析方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）極差分析：就是將每列不同水平對應(yīng)的試驗(yàn)結(jié)果進(jìn)行求和，然后對求和結(jié)果計(jì)算極差，并找出最大值所在的水平；極差最大所在因素,以此表示因素不同水平的差異性。計(jì)算發(fā)現(xiàn)：因素A的極差最大因素C次之因素B最小結(jié)論：因素的主次順序

ACB問題：它們的影響是否顯著呢？6.4.2正交試驗(yàn)例子分析—方差分析方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）

自由度

rt-1=n-1

r-1

n-1-m(r-1)說明：m*為非空列個數(shù)(即實(shí)際因素個數(shù)),如果因素與誤差列合并，則其自由度相應(yīng)增加。6.4.2正交試驗(yàn)例子分析—方差分析方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）說明:此處空白列作為誤差E，又因?yàn)橐蛩谺的平方和小于誤差E,所以將因素B和誤差E合并作為誤差E，然后建立F統(tǒng)計(jì)量用來檢驗(yàn)因素的顯著性。注意：此時誤差列的自由度由2變?yōu)?。問題：如果沒有空白列做誤差，怎么分析？6.4.2正交試驗(yàn)例子分析—R過程R中的包：DoE.base,提供函數(shù)oa.design可生成正交表。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）oa.design(ID=NULL,nruns=NULL,nfactors=NULL,nlevels=NULL,s,columns="order",repeat.only=FALSE,randomize=TRUE,seed=NULL,…)其中ID是試驗(yàn)表格的編號，如L9.3.4，L8.2.7等，可以通過show.oas()函數(shù)得到；nruns表示試驗(yàn)次數(shù)，nfactors表示因素個數(shù)，nlevels表示水平個數(shù)，s表示因素的名稱，randomize表示試驗(yàn)順序是否隨機(jī)產(chǎn)生，seed表示隨機(jī)數(shù)的種子.LT=oa.design(nfactors=4,nlevels=3,s=c('溫度','時間','紙漿濃度','空列'),randomize=FALSE)6.4.2正交試驗(yàn)例子分析—R過程方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）show.oas(name="all",nruns="all",nlevels="all",factors="all",show=10,parents.only=FALSE)

主要用來查詢正交表show.oas(show=20)顯示前20個正交表show.oas(nlevels=c(2,3)顯示水平含為2,3的正交表show.oas(factors=list(nlevels=c(2,3),number=c(1,3)))

顯示支持1個2水平，3個3水平的正交表4.show.oas(nruns=c(4,20))顯示試驗(yàn)次數(shù)介于4到20的正交表6.4.2正交試驗(yàn)例子分析—R過程方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）install.packages(‘DoE.base’);library(DoE.base)LT=oa.design(nfactors=4,nlevels=3,s=c('溫度','時間','紙漿濃度','空列'),randomize=F)responses=c(3150,3030,3100,2830,3160,2950,2910,2520,2670)#變量名必須取名為responsesL9=add.response(LT,responses)#加上一個響應(yīng)變量,完整的正交+試驗(yàn)表aov.L9=aov(responses~溫度+時間+紙漿濃度,data=L9);summary(aov.L9)#發(fā)現(xiàn)時間因素的均方和小于空列，故將時間因素與空列合并成誤差aov.L9=aov(responses~溫度+紙漿濃度,data=L9);summary(aov.L9)6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）

有交互作用的試驗(yàn)，情況比較復(fù)雜，選擇正交表和安排因素所在列都需要認(rèn)真考慮。所謂因素之間的交互作用指的是因素之間的不同水平搭配對試驗(yàn)結(jié)果的影響。我們只考慮二級交互作用，即二因素之間的交互作用。例如:A:B,A:C,B:C分別表示A和B之間,A和C之間,B和C之間的交互作用。而A,B,C各自的作用稱為主效應(yīng)。6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）

只有配備交互作用列表的正交表才能安排交互作用，才能確定任意兩個因素交互作用所在的列，這種安排過程稱為表頭設(shè)計(jì)。

表中的元素為交互作用所在的列號，

比如第1列和第2列的交互作用列，應(yīng)查表中第1行第2列的元素：3，它表示該交互作用列應(yīng)安排在第3列；

同理第3列與第5列的交互作用，即查表中第3行第5列的元素：6，表示該交互作用列應(yīng)安排在第6列，其余類推。

6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）說明:此處空白列作為誤差E,交互作用列保留數(shù)字即可。說明:直觀分析同無交互模型6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）說明:優(yōu)選方案A2(真優(yōu)11)B2(8萬苗/畝)C2(15公斤)--極差分析

因子主次：BAXBACBXCAXC6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）因子主次：BAXBA方法1全面搭配比較方法2最優(yōu)搭配比較6.4.3正交試驗(yàn)—有交互作用方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）

從方差分析表中，我們看到因素B和交互作用AXB對指標(biāo)的影響極其顯著；因素A，C比較顯著；而交互AXC和BXC對指標(biāo)的影響不顯著。此處：剩余誤差的自由度為2，為什么？6.2協(xié)方差分析(補(bǔ)充-不考核)為了提高試驗(yàn)的精確性和準(zhǔn)確性，我們對除研究因素以外的一切條件都需要采取有效措施嚴(yán)加控制，使它們在因素的不同水平間盡量保持一致，這叫做試驗(yàn)控制。但當(dāng)我們進(jìn)行試驗(yàn)設(shè)計(jì)時，即使做出很大努力控制，也經(jīng)常會碰到試驗(yàn)個體的初始條件不同的情況，如果不考慮這些因素有可能導(dǎo)致結(jié)果失真。如果考慮這些不可控的因素，這種方差分析就叫做協(xié)方差分析，其是將回歸分析和方差分析結(jié)合在一起的方法。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2協(xié)方差分析(補(bǔ)充-不考核)基本原理如下:將一些對響應(yīng)變量Y有影響的變量X(未知或難以控制的因素)看作協(xié)變量，建立響應(yīng)變量Y隨X變化的線性回歸分析，從Y的總的平方和中扣除X對Y的回歸平方和，對殘差平方和作進(jìn)一、步分解后再進(jìn)行方差分析。比如同一因素不同水平起點(diǎn)不一的前提下所做的實(shí)驗(yàn)的分析因素的不同水平作為協(xié)變量。有些情況下，初始狀態(tài)時各水平的值不在同一起點(diǎn)，這是人為無法控制的，因此將此也作為一個協(xié)變量。方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2協(xié)方差分析(補(bǔ)充)install.packages(‘HH’)library(HH)Weight_Initial=c(15,13,11,12,12,16,14,17,17,16,18,18,21,22,19,18,22,24,20,23,25,27,30,32)Weight_Increment=c(85,83,65,76,80,91,84,90,97,90,100,95,103,106,99,94,89,91,83,95,100,102,105,110)feed=gl(3,8,24)data_feed=data.frame(Weight_Initial,Weight_Increment,feed)summary(ancova(Weight_Increment~Weight_Initial+feed,data=data_feed))Summary(ancova(Weight_Increment~Weight_Initial*feed,data=data_feed))方法篇：方差分析與正交試驗(yàn)設(shè)計(jì)(第六章）6.2協(xié)方差分析(補(bǔ)充)DfSumSqMeanSqFvaluePr(>F)Weight_Initial11621.11621.1142.441.50e-10***feed2707.2353.631.077.32e-07***Residuals20227.611.4--