R-多元統(tǒng)計分析上機講義

應(yīng)用多元統(tǒng)計分析AppliedMultivariateStatisticalAnalysis第一章緒論在實際問題中，很多隨機現(xiàn)象涉及到的變量不是一個，而是經(jīng)常是多個變量，并且這些變量間又存在一定的聯(lián)系。我們經(jīng)常需要處理多個變量的觀測數(shù)據(jù)，如果用一元統(tǒng)計方法由于忽視了各個變量之間可能存在的相關(guān)性，一般說來，丟失信息太多，分析的結(jié)果不能客觀全面反映數(shù)據(jù)所包含的內(nèi)容，因此，我們就需要用到多元統(tǒng)計的方法。多元統(tǒng)計分析(MultivariateStatisticalAnalysis)也稱多變量統(tǒng)計分析、多因素統(tǒng)計分析或多元分析，是研究客觀事物中多變量(多因素或多指標)之間的相互關(guān)系和多樣品對象之間差異以及以多個變量為代表的多元隨機變量之間的依賴和差異的現(xiàn)代統(tǒng)計分析理論和方法。多元統(tǒng)計分析是解決實際問題的有效的數(shù)據(jù)處理方法。隨著電子計算機使用的日益普及，多元統(tǒng)計統(tǒng)計方法已廣泛地應(yīng)用于自然科學(xué)、社會科學(xué)的各個方面。第二章矩陣矩陣即是二維的數(shù)組，它非常的重要，以至于需要單獨討論。由于矩陣應(yīng)用非常廣泛，因此對它定義了一些特殊的應(yīng)用和操作，R包括許多只對矩陣操作的操作符和函數(shù)。矩陣的建立在R中最為常用的是用命令matrix()建立矩陣，而對角矩陣常用函數(shù)diag()建立例如>X<-matrix(1,nr=2,nc=2)>X[,1][,2]TOC\o"1-5"\h\z[1,]11[2,]11X<-diag(3)#生成單位陣X[,1][,2][,3][1,]100[2,]010[3,]001>diag(2.5,nr=3,nc=5)[,1][,2][,3][,4][,5][1,]2.50.00.000[2,]0.02.50.000[3,]0.00.02.500>X<-matrix(1:4,2)#等價于X<-matrix(1:4,2,2)>X[,1][,2][1,]13[2,]24>rownames(X)<-c("a","b")>colnames(X)<-c("c","d")>Xcda13b24>dim(X)[1]22>dimnames(X)[[1]][1]"a""b"[[2]][1]"c""d"注意：循環(huán)準則仍然適用于matrix（），但要求數(shù)據(jù)項的個數(shù)等于矩陣的列數(shù)的倍數(shù)，否則會出現(xiàn)警告。矩陣的維數(shù)使用c（）會得到不同的結(jié)果（除非是方陣），因此需要小心。數(shù)據(jù)項填充矩陣的方向可通過參數(shù)byrow來指定，其缺省是按列填充的（byrow二FALSE），byrow二TRUE表示按行填充數(shù)據(jù)。再看幾個例子：X<-matrix（1:4，2，4）#按列填充TOC\o"1-5"\h\zX［，1］［，2］［，3］［，4］［1，］1313［2，］2424X<-matrix（1:4，2，3）Warningmessage:Inmatrix（1:4，2，3）:數(shù)據(jù)長度［4］不是矩陣列數(shù)［3］的整倍數(shù)X<-matrix（1:4，c（2，3））#不經(jīng)常使用X［，1］［，2］［1，］13［2，］24X<-matrix（1:4，2，4，byrow=TRUE）#按行填充X[,1][,2][,3][,4][1,]1234[2,]1234因為矩陣是數(shù)組的特例,R中數(shù)組由函數(shù)array()建立，因此矩陣也可以用函數(shù)array()來建立，其一般格式為：array(data,dim,dimnames)其中data為一向量，其元素用于構(gòu)建數(shù)組；dim為數(shù)組的維數(shù)向量(為數(shù)值型向量)；dimnames為由各維的名稱構(gòu)成的向量(為字符型向量)，缺省為空。看幾個例子：A<-array(1:6,c(2,3))A[,1][,2][,3][1,]135[2,]246A<-array(1:4,c(2,3))A[,1][,2][,3][1,]131[2,]242A<-array(1:8,c(2,3))A[,1][,2][,3][1,]135[2,]2462.2矩陣的下標(index)與子集(元素)的提取矩陣的下標可以使用正整數(shù)、負整數(shù)和邏輯表達式，從而實現(xiàn)子集的提取或修改。考查矩陣x<-matrix(1:6,2,3)x[,1][,2][,3][1,]135[2,]246?提取一個元素>x[2,2][1]4?提取若一個或若干個行或列>x[2,2][1]4>x[2,][1]246>x[,2][1]34>x[,2,drop=FALSE][,1][1,]3[2,]4>x[,c(2,3),drop=FALSE][,1][,2][1,]35[2,]46?去掉某一個或若干個行與列>x[-1,][1]246>x[,-2][,1][,2][1,]15[2,]26?添加與替換元素>x[,3]<-NA>x[,1][,2][,3][1,]13NA[2,]24NA>x[is.na(x)]<-1#缺失值用1代替>x[,1][,2][,3][1,]131[2,]241矩陣四則運算矩陣也可以進行四則運算（“+”、“-”、“*”、“/”，“’”），分別解釋為矩陣對應(yīng)元素的四則運算。在實際應(yīng)用中，比較有實際應(yīng)用的是矩陣的相加，相減，相乘和矩陣的求逆。矩陣的加減運算一般要求矩陣形狀完全相同（dim屬性完全相同），矩陣的相乘一般要求一矩陣的列維數(shù)與另一矩陣的行維數(shù)相同，而矩陣要求逆的話，一般要求它為一方陣。矩陣的加減運算若A，B為兩個形狀相同的矩陣，兩矩陣的和為C，R中表達式為：C<-A+B兩矩陣的差為D，R中表達式為：D<-A-B矩陣也可以與數(shù)進行加減，A+5表示A中的每個元素加上5。矩陣的相乘操作符％*%用于矩陣相乘。若矩陣A的列數(shù)等于矩陣B的行數(shù)，矩陣A乘以矩陣B表示為：A%*%B注：X*Y表示兩個矩陣的逐元相乘，而不是X和Y的乘積。矩陣的求逆若矩陣A為一方陣，矩陣的逆可以用下面的命令計算：solve(A)。操作符solve()可以用來求解線性方程組：Ax=b，解為solve(A,b)在數(shù)學(xué)上，用直接求逆的辦法解x<-solve(A)%*%b相比solve(A,b)不僅低效而且還有一種潛在的不穩(wěn)定性。矩陣的其他一些代數(shù)運算求轉(zhuǎn)置矩陣轉(zhuǎn)置函數(shù)為t()，矩陣X的轉(zhuǎn)置為t(X)。提取對角元素提取對角元的函數(shù)為diag()。例如：X<-matrix(1:4,2,2)diag(X)[1]14事實上，diag()的作用依賴于自變量，diag(vector)返回以自變量(向量)為主對角元素的對角矩陣；diag(matrix)返回由矩陣的主對角元素所組成的向量；diag(k)(k為標量)返回k階單位陣。矩陣的合并與拉直函數(shù)cbind()把幾個矩陣橫向拼成一個大矩陣，這些矩陣行數(shù)應(yīng)該相同；函數(shù)rbind()把幾個矩陣列向拼成一個大矩陣，這些矩陣列數(shù)應(yīng)該相同。(如果參與合并的矩陣比其它矩陣行數(shù)少或列數(shù)少，則循環(huán)不足后合并。)例如：m1<-matrix(1,nr=2,nc=2)TOC\o"1-5"\h\zm1[,1][,2][1,]11[2,]11>m2<-matrix(2,nr=2,nc=2)>m2[,1][,2]TOC\o"1-5"\h\z[1,]22[2,]22>rbind(m1,m2)[,1][,2][1,]11[2,]11[3,]22[4,]22>cbind(m1,m2)[,1][,2][,3][,4][1,]1122[2,]1122方陣的行列式求方陣的行列式使用det(：X<-matrix(1:4,2)>X[,1][,2][1,]13[2,]24>det(X)[1]-2矩陣的特征根和特征向量函數(shù)eigen()用來計算矩陣的特征值和特征向量。這個函數(shù)的返回值是一個含有values和vectors兩個分量的列表。命令A(yù)<-eigen(X)>A$values[1]5.3722813-0.3722813$vectors[,1][,2][1,]-0.5657675-0.9093767[2,]-0.82456480.4159736MatrixfacilitesInthefollowingexamples,AandBarematricesandxandbareavectors.OperatororFunctionDescriptionA*BElement-wisemultiplicationA%*%BMatrixmultiplicationA%o%BOuterproduct.AB'crossprod(A,B)crossprod(A)A'BandA'Arespectively.t(A)Transposediag(x)Createsdiagonalmatrixwithelementsofxintheprincipaldiagonaldiag(A)Returnsavectorcontainingtheelementsoftheprincipaldiagonaldiag(k)Ifkisascalar,thiscreatesakxkidentitymatrix.Gofigure.solve(A,b)Returnsvectorxintheequationb=Ax(i.e.,A-1b)solve(A)InverseofAwhereAisasquarematrix.ginv(A)Moore-PenroseGeneralizedInverseofA.ginv(A)requiresloadingtheMASSpackage.y<-eigen(A)y$valaretheeigenvaluesofAy$vecaretheeigenvectorsofAy<-svd(A)SinglevaluedecompositionofA.y$d=vectorcontainingthesingularvaluesofAy$u=matrixwithcolumnscontaintheleftsingularvectorsofAy$v=matrixwithcolumnscontaintherightsingularvectorsofAR<-chol(A)CholeskifactorizationofA.Returnstheuppertriangularfactor,suchthatR'R=A.y<-qr(A)QRdecompositionofA.y$qrhasanuppertrianglethatcontainsthedecompositionandalowertrianglethatcontainsinformationontheQdecomposition.y$rankistherankofA.y$qrauxavectorwhichcontainsadditionalinformationonQ.y$pivotcontainsinformationonthepivotingstrategyused.cbind(A,B,...)Combinematrices(vectors)horizontally.Returnsamatrix.rbind(A,B,...)Combinematrices(vectors)vertically.Returnsamatrix.rowMeans(A)Returnsvectorofrowmeans.rowSums(A)Returnsvectorofrowsums.colMeans(A)Returnsvectorofcolumnmeans.colSums(A)Returnsvectorofcoumnsums.2.4.7.其它函數(shù)交叉乘積(crossproduct),函數(shù)為crossprod(),crossprod(X,Y)表示一般的內(nèi)積XzY,即X的每一列與Y的每一列的內(nèi)積組成的矩陣；QR分解，函數(shù)為qr()，矩陣X的QR分解為X=QR,Q為正交陣，R為上三角陣；等等。矩陣的統(tǒng)計運算函數(shù)cov()和cor()分別用于計算矩陣的協(xié)方差陣和相關(guān)系數(shù)陣。矩陣的排列是有方向性的，在R中規(guī)定矩陣是按列排的，若沒有特別說明，函數(shù)max()，min()，median()，var()，sd(),sum(),cumsum(),cumprod(),cummax(),cummin()的使用對于矩陣也是按列計算的，但也可以通過選項MARGIN來改變。下面我們要用到對一個對象施加某種運算的函數(shù)apply()，其格式為apply(X,MARGIN,FUN)其中X為參與運算的矩陣,FUN為上面的一個函數(shù)或“+”、“_”、“*”、“\”(必須放在引號中)，MARGIN=1表示按列計算，MARGIN=2表示按行計算。我們還用到sweep()函數(shù)，命令sweep(X,MARGIN,STATS,FUN)表示從矩陣X中按MATGIN計算STATS，并從X中除去(sweepout)。求均值m<-matrix(rnorm(n=12),nrow=3)apply(m,MARGIN=1,FUN=mean)#求各行的均值-0.37738650.38641380.2052353apply(m,MARGIN=2,FUN=mean)#求各列的均值[1]0.33862020.7320669-0.4624578-0.3225460標準化scale(m,center=T,scale=T)減去中位數(shù)row.med<-apply(m,MARGIN=1,FUN=median)>sweep(m,MARGIN=1,STATS=row.med,FUN=”-”)第三章多元正態(tài)分布及參數(shù)的估計繪制二元正態(tài)密度函數(shù)及其相應(yīng)等高線圖書上例2.2.，2錯誤！未找到引用源。時的二元正態(tài)密度函數(shù)及其等高線圖：x<-seq(-3,3,by=0.1)y<-xf<-function(x,y,a=1,b=1,r=0){a1=sqrt(a)b1=sqrt(b)d=1-r*rd1=sqrt(d)*a1*b1z=1/(2*pi*d1)*exp((-x*x/a-y*y/b+2*r*x*y/(a1*b1))/(2*d))}z<-outer(x,y,f)#外積函數(shù)persp(x,y,z,xlim=range(x),ylim=range(y),zlim=range(z,na.rm=TRUE),theta=30,nticks=5,ticktype="detailed",sub="O1=02=1,p=0時的二元正態(tài)密度函數(shù)")#密度函數(shù)圖contour(x,y,z)#等高線圖image(x,y,z)#等高線圖,實際數(shù)據(jù)大小用不同色彩表示所得圖形為：

cf1=ct2=1,卩二0時的二元正態(tài)密度函數(shù)相應(yīng)等高線圖Outer(x,y,f)是一個一般性的外積函數(shù)，調(diào)用函數(shù)f,把x的任一個元素與y的任意一個元素搭配起來作為f的自變量計算得到新的元素值，當函數(shù)缺省時表示乘積情況。對參數(shù)進行修改，可以繪制任一二元正態(tài)密度函數(shù)及其相應(yīng)的等高線圖。多元正態(tài)分布的參數(shù)估計多元正態(tài)總體的相關(guān)量設(shè)觀測數(shù)據(jù)陣為?樣本均值向量設(shè)錯誤！未找到引用源。，錯誤！未找到引用源。=1，2，錯誤！未找到引用源。，錯誤！未找到引用源。，則樣本均值向量Xn：錯誤！未找到引用源。，由2.5.1可得：Xn<-apply(x,MARGIN=2,mean)或者ln<-rep(1，n)Xn<-(ln%*%x)/nXn即為所求樣本均值向量。樣本離差陣(交叉乘積陣)樣本離差陣A：錯誤！未找到引用源。。A<-crossprod(x)-2*Xn%*%t(Xn)或者m<-diag(1，n)-matrix(1，n，n)/nA<-t(x)%*%m%*%xA即為所求樣本離差陣。樣本協(xié)方差陣R中求樣本協(xié)方差陣的函數(shù)為cov()。樣本數(shù)據(jù)陣X的協(xié)方差矩陣S即為:S<-cov(X)?樣本相關(guān)陣R中求樣本協(xié)方差陣的函數(shù)為cor()。樣本數(shù)據(jù)陣X的協(xié)方差矩陣R即為:>R<-cor(X)極大似然估計極大似然估計法是建立在極大似然原理基礎(chǔ)上的一種統(tǒng)計方法。設(shè)總體X，其概率密度函數(shù)(連續(xù)情況)或分布律(離散情況)為錯誤！未找到引用源。，其中錯誤！未找到引用源。是未知參數(shù)(或未知參數(shù)向量)。設(shè)X],X2,…，Xn為取自總體X的樣本，貝y似然函數(shù)錯誤！未找到引用源。為：錯誤！未找到引用源?！?錯誤！未找到引用源。)=錯誤！未找到引用源。求使似然函數(shù)達到最大的參數(shù)錯誤！未找到引用源。的值，即極大似然估計值。在單參數(shù)場合，在R中可以使用函數(shù)optimize()求極大似然估計值。optimize()的調(diào)用格式如下：optimize(f=,interval=,lower=min(interval),

upper=max(interval),maximum=TRUE,

tol=.Machine$double.eps"0.25,…)說明：f是似然函數(shù)，interva1是參數(shù)錯誤！未找到引用源。的取值范圍，lower是錯誤!未找到引用源。的下界，upper是錯誤！未找到引用源。的上界，maximum=TRUE是求極大值，否則(maximum=FALSE)表示求函數(shù)的極小值，tol是表示求值的精確度，…是對f的附加說明。在多參數(shù)場合，在R中用函數(shù)optim(咸者nlm()來求似然函數(shù)的極大值，并求相應(yīng)的極大值點。optim()的調(diào)用格式如下：optim(par,fn,gr=NULL,method=c("Nelder-Mead","BFGS","CG","L-BFGS-B","SANN"),lower=-Inf,upper=Inf,control=list(),hessian=FALSE,…)nlm()的定義如下：nlm(f,p,hessian=FALSE,typsize=rep(1,length(p)),fscale=1,

print.level=0,ndigit=12,gradtol=1e-6,stepmax=max(1000*sqrt(sum((p/typsize)"2)),1000),steptol=1e-6,iterlim=100,check.analyticals=TRUE,…)三者主要區(qū)別是：函數(shù)nlm()僅使用牛頓-拉夫遜算法求函數(shù)的最小值點；函數(shù)optim()提供method選項給出的5種方法中的一種進行優(yōu)化；上面二個可用于多維函數(shù)的極值問題,而函數(shù)optimize()僅適用于一維函數(shù)，但可以用于最大與最小值點。(具體選項見幫助。)第四章多元正態(tài)總體參數(shù)的假設(shè)檢驗在一元統(tǒng)計中，用于檢驗一元正態(tài)總體參數(shù)錯誤！未找到引用源。，錯誤！未找到引用源。的抽樣分布有錯誤！未找到引用源。分布，錯誤！未找到引用源。分布、F分布風(fēng)，它們都是來自總體錯誤！未找到引用源。的隨機樣本導(dǎo)出的檢驗統(tǒng)計量。推廣到多元正態(tài)總體后，也有相應(yīng)于以上三個常用分布的統(tǒng)計量：威沙特(Wishart)統(tǒng)計量，霍特林(Hotelling)錯誤！未找到引用源。統(tǒng)計量，威爾克斯(Wilks)錯誤！未找到引用源。統(tǒng)計量，這些統(tǒng)計量是多元統(tǒng)計分析所涉及的假設(shè)檢驗問題的基礎(chǔ)。幾個重要統(tǒng)計量的分布對于多元正態(tài)總體來說，存在幾個重要的統(tǒng)計量：威沙特(Wishart)統(tǒng)計量，霍特林(Hotelling)錯誤！未找到引用源。統(tǒng)計量，威爾克斯(Wilks)錯誤！未找到引用源。統(tǒng)計量等,討論這些統(tǒng)計量的分布是多元統(tǒng)計分析所涉及的假設(shè)檢驗問題的基礎(chǔ)。單總體均值向量的檢驗及置信域均值向量的檢驗書上例3.2.1，R程序如下>x<-matrix(c(3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.1,7.6,6.7,47.4,8.5,5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,5.5,40.9,9.4),20,3,byrow=TRUE)>n<-20>p<-3>u0<-c(4,50,10)#所給總體均值>ln<-rep(1,20)>x0<-(ln%*%x)/n#樣本均值>xm<-x0-u0>mm<-diag(1,20)-matrix(1,20,20)/n>a<-t(x)%*%mm%*%x#樣本離差陣>ai=solve(a)>dd=xm%*%ai%*%t(xm)>d2=(n-1)*dd>t2=n*d2;f<-(n-p)*t2/((n-1)*p)#檢驗統(tǒng)計量f[,1][1,]2.904546fa〈-qf(0.95,p,n-p)#自由度為(p,n-p)的F分布的0.95分位數(shù)fa[1]3.196777b〈-1-pf(f,p,n-p)#尾概率值b[,1][1,]0.06492834beta〈-pf(fa,p,n-p,t2)#犯第二類錯誤的概率(假設(shè)總體均值錯誤！未找到引用源。)beta[1]0.3616381取檢驗水平為錯誤！未找到引用源。0.05,由尾概率值p=0.06492834錯誤！未找到引用源。0.05=錯誤！未找到引用源。，可得錯誤！未找到引用源。相容；同樣由F=2.904546錯誤！未找到引用源。3.196777=Fa，也可得錯誤！未找到引用源。相容。在這種情況下，可能犯第二類錯誤,概率為錯誤！未找到引用源。=0.3616(假定總體均值錯誤！未找到引用源。)。樣本協(xié)方差陣的特征值和特征向量書上例3.2.2，R程序為：x〈-matrix(c(3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.1,7.6,6.7,47.4,8.5,5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,5.5,40.9,9.4),20,3,byrow=TRUE)>s<-cov(x)>s[,1][,2][,3][1,]2.87936810.0100-1.809053[2,]10.010000199.7884-5.640000[3,]-1.809053-5.64003.627658>a<-eigen(s)$values[1]200.4624644.5315911.301392$vectors[,3][,1][,2][1,]-0.05084144-0.573703640.81748351[2,]-0.998283520.05302042-0.02487655[3,]0.029071560.817345080.57541452多總體均值向量的檢驗兩正態(tài)總體均值向量的檢驗書上例3.3.1,R程序為：>n<-10>m<-10>p<-4>x<-matrix(c(65,75,60,75,70,55,60,65,60,55,+35,50,45,40,30,40,45,40,50,55,25,20,35,40,30,+35,30,25,30,35,60,55,65,70,50,65,60,60,70,75),10)>ln<-rep(1,n)>x0<-(ln%*%x)/n>mx<-diag(1,n)-matrix(1,n,n)/n>a1<-t(x)%*%mx%*%x>y<-matrix(c(55,50,45,50,55,60,65,50,40,45,+55,60,45,50,50,40,55,60,45,50,40,45,35,50,+30,45,45,35,30,45,65,70,75,70,75,60,75,80,65,70),10)>y0<-(ln%*%y)/n>my<-diag(1,n)-matrix(1,n,n)/n>a2<-t(y)%*%my%*%y>a<-a1+a2>xy<-x0-y0>ai<-solve(a)

>dd<-xy%*%ai%*%t(xy)>d2<-(m+n-2)*dd>t2<-n*m*d2/(n+m)>f<-(n+m-1-p)*t2/((n+m-2)*p)>pp<-1-pf(f,p,m+n-p-1)>x0[,1][,2][,3][,4][1,]644330.563>y0[,1][,2][,3][,4][1,]51.5514070.5>a1[,1][,2][,3][,4][1,]490-170-120.0-245[2,]-17051010.0310[3,]-12010322.5260[4,]-245310260.0510>a2[,1][,2][,3][,4][1,]502.560175-7.5[2,]60.039050195.0[3,]175.050450-100.0[4,]-7.5195-100322.5>d2[,1][1,]5.972499>t2[,1][1,]29.86250>f[,1][1,]6.221353>pp[,1][1,]0.003705807取檢驗水平為錯誤！未找到引用源。0.01，根據(jù)尾概率值p=0.003705807錯誤！未找到引用源。0.01=錯誤！未找到引用源。，可得應(yīng)否定錯誤！未找到引用源。。多個正態(tài)總體均值向量的檢驗-多元方差分析書上例3.3.2，可利用類似例3.2.1或例3.3.1的程序進行計算得出結(jié)論。下面我們用R自帶的manova()函數(shù)進行分析。程序如下：

x<-read.table("D:/data/d332.txt",header=T)x<-as.matrix(x[,1:4])rate<-factor(gl(3,20),labels=c("group1","group2","group3"))fit<-manova(x~rate)summary.aov(fit)#對每一個變量進行單因素方差分析summary(fit,test="Wilks")#使用威爾克斯錯誤！未找到引用源。統(tǒng)計量程序結(jié)果：>summary.aov(fit)Responsex1:SumSqMeanSqFvaluePr(>F)0.0004401***rateResidualsDf257390661254091953322008.878Signif.codes:0‘***'0.001‘**'0.01‘*'0.05‘.Responsex2:DfSumSqMeanSqFvaluePr(>F)rate2401720092.82930.06738.Residuals5740467710Signif.codes:0‘***'0.001‘**'0.01‘*'0.05‘.Responsex3:DfSumSqMeanSqFvaluePr(>F)rate213.436.720.18380.8326Residuals572082.5036.54Responsex4:DfSumSqMeanSqFvaluePr(>F)rate217.208.600.47850.6222Residuals571024.4017.97>summary(fit,test="Wilks")DfWilksapproxFnumDfdenDfPr(>F)rate20.662123.0906981080.003538**Residuals57Signif.codes:0‘***'0.001‘**'0.01‘*'0.05‘.0.10.10.1結(jié)果說明：(1)取檢驗水平為0.01,則對四個指標逐項用一元方差分析方法進行檢驗，由p值可得三個組指標間只有第一個指標錯誤！未找到引用源。有顯著差異(錯誤！未找到引用源。=0.0004401)；(2)取檢驗水平為0.01，利用威爾克斯錯誤！未找到引用源。統(tǒng)計量得到p=0.003538錯誤！未找到引用源。0.01,故拒絕原假設(shè)，即認為三個組的指標之間有顯著差異。協(xié)方差陣的檢驗多總體協(xié)方差陣的檢驗書上例3.4.1，R程序略(類似例3.2.1或例3.3.1)獨立性檢驗書中例3.5.1，R程序為：x<-matrix(c(3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,+3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.1,7.6,6.7,47.4,8.5,+5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,+1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,+5.5,40.9,9.4),20,3,byrow=TRUE)n<-20p<-3x0<-(ln%*%x)/n>xm<-x0-u0>mm<-diag(1,20)-matrix(1,20,20)/n>a<-t(x)%*%mm%*%x>a0<-det(a)a1<-a[1,1]*a[2,2]*a[3,3]v<-a0/a1b<-n-1.5-(p*p*p-3)/(3*p*p-3*3)df<-0.5*(p*(p+1)-2*3)kc<--b*log(v)p0<-1-pchisq(kc,df)kc[1]9.755514>p0[1]0.02076288取檢驗水平為錯誤！未找到引用源。0.05,根據(jù)尾概率值p=0.02076288錯誤！未找到引用源。0.05=錯誤！未找到引用源。，可得應(yīng)否定錯誤！未找到引用源。，由R軟件所的結(jié)果與SAS軟件所的結(jié)果一致。正態(tài)性檢驗書中例3.6.1，R程序為：>x<-matrix(c(100,99,96,99,96,75,97,68,76,62,67,34,100,97,100,96,78,97,+89,88,84,39,78,37),12)>n<-12>p<-2>ln<-rep(1,n)>x0<-(ln%*%x)/n>s<-cov(x)>si<-solve(s)>m<-0>for(iin1:n){+