R語(yǔ)言之方差分析篇

R語(yǔ)?之?差分析篇當(dāng)包含的因?是解釋變量時(shí)，通常會(huì)從預(yù)測(cè)轉(zhuǎn)向級(jí)別差異的分析，即稱(chēng)作?差分析（ANOVA）組間因?因變量?變量均衡設(shè)計(jì)（balanceddesign）組內(nèi)因?單因素組間?差分析單因素組內(nèi)?差分析重復(fù)測(cè)量?差分析主效應(yīng)交叉效應(yīng)因素?差分析混淆因素?擾變數(shù)協(xié)變量協(xié)?差1、ANOVA模型擬合1.aov（）函數(shù)語(yǔ)法：aov（formula,data=dataframe）R表達(dá)式中的特殊符號(hào)符號(hào)?法~分隔符號(hào)，左邊為響應(yīng)變量，右邊為解釋變量eg：y~A+B+C+分隔解釋變量：表?變量的交互項(xiàng)eg：y~A+B+A:B*表?所有可能交互項(xiàng)eg：y~A*B*C可展開(kāi)為：y~A+B+C+A:B+A:C+B:C+A:B:C^表?交互項(xiàng)達(dá)到次數(shù)eg：y~(A+B+C)^2展開(kāi)為：y~A+B+C+A:B+A:C+B:C.表?包含除因變量外的所有變量eg：若?個(gè)數(shù)據(jù)框包括變量y,A、B和C，代碼y~.可展開(kāi)為y~A+B+C常見(jiàn)研究設(shè)計(jì)的表達(dá)式設(shè)計(jì)表達(dá)式單因素ANOVAy~A含單個(gè)協(xié)變量的單因素ANCOVAy~x+A雙因素ANOVAy~A*B含兩個(gè)協(xié)變量的雙因素ANCOVAy~x1+x2+A*B隨機(jī)化區(qū)組y~B+A（B是區(qū)組因?）單因素組內(nèi)ANOVAy~A+Error（Subject/A）含單個(gè)組內(nèi)因?（W）和單個(gè)組間因?（B）的重復(fù)測(cè)量ANOVAy~B*W+Error（Subject/W）?平衡設(shè)計(jì)時(shí)或存在協(xié)變量時(shí)，效應(yīng)項(xiàng)的順序?qū)Y(jié)果影響較?越基礎(chǔ)的效應(yīng)應(yīng)越需要放在表達(dá)式前?，?先是協(xié)變量、然后是主效應(yīng)、接著是雙因素的交互項(xiàng)，再接著是三因素的交互項(xiàng)若研究不是正交的，?定要謹(jǐn)慎設(shè)置疚的順序2、單因素?差分析（1）單因素?差分析#單因素?差分析（感興趣地是?較分類(lèi)因?定義的兩個(gè)或多個(gè)組別中的因變量均值）install.packages("multcomp")library(multcomp)attach(cholesterol)str(cholesterol)cholesteroltable(trt)aggregate(response,by=list(trt),FUN=mean)aggregate(response,by=list(trt),FUN=sd)fit<-aov(response~trt)summary(fit)library(gplots)plotmeans(response~trt,xlab="Treatment",ylab="Response",main="MeanPlot\nwith95%CI")（2）多重?較多重?較?于解決某?組別與其他的不同TukeyHSD（）函數(shù)提供了對(duì)各組均值差異的成對(duì)檢驗(yàn)，但與HH包存在兼容性問(wèn)題(（某些版本中）；TukeyHSD(fit)par(las=2)par(mar=c(5,4,6,2))plot(TukeyHSD(fit))library(multcomp)par(mar=c(5,4,6,2))tuk<-glht(fit,linfct=mcp(trt="Tukey"))plot(cld(tuk,level=0.05),col="lightgrey")（3）評(píng)估檢驗(yàn)的假設(shè)條件當(dāng)因變量服從正態(tài)頒，各組?差相等時(shí)，可?Q-Q圖來(lái)檢驗(yàn)正態(tài)性假設(shè)qqPlot（）要求?lm（）擬合，若數(shù)據(jù)落在95%的置信區(qū)間范圍內(nèi)，說(shuō)明滿?正態(tài)性假設(shè)。R提供的可以做?差齊性檢驗(yàn)的函數(shù)Bartlett檢驗(yàn)bartlet.test（）Fligner-Killeen檢驗(yàn)fligner.test（）Brown-Forsythe檢驗(yàn)離群點(diǎn)檢驗(yàn)car包中的outlierTest（）函數(shù)來(lái)檢測(cè)離群點(diǎn)3、單因素協(xié)?差分析install.packages("multcomp")library(multcomp)head(litter,n=21)data(litter,package="multcomp")attach(litter)options(digits=5)table(litter$dose)aggregate(weight,by=list(litter$dose),FUN=mean)fit<-aov(weight~gesttime+litter$dose)summary(fit)因使?了協(xié)變量，短途運(yùn)輸獲取調(diào)整的組均值即去除協(xié)變量疚后的組均值，可使?effects包中的effects（）函數(shù)來(lái)計(jì)算調(diào)整的均值library(effects)effect("dose",fit)?戶定義的對(duì)照的多重?較library(multcomp)contrast<-rbind("nodrugvs.drug"=c(3,-1,-1,-1))summary(glht(fit,linfct=mcp(dose=contrast)))（1）評(píng)估檢驗(yàn)的假設(shè)條件ANCOVA與ANOVA相同，都城要正態(tài)性和同?差性假設(shè)另ANOCVA還假定回歸低低斜率相同，eg當(dāng)ANCOVA模型飲?懷孕時(shí)間*劑量的交互項(xiàng)時(shí)，可對(duì)回歸斜率的同質(zhì)性進(jìn)?檢驗(yàn)。eg檢驗(yàn)回歸斜率的同質(zhì)性library(multcomp)fit2<-aov(weight~gesttime*dose,data=litter)summary(fit2)（2）結(jié)果的可視化install.packages("HH")library(HH)ancova(weight~gesttime+dose,data=litter)4、雙因素?差分析<strong>雙因素ANOVA</strong>attach(ToothGrowth)head(ToothGrowth)table(supp,ToothGrowth$dose)aggregate(len,by=list(supp,ToothGrowth$dose),FUN=mean)aggregate(len,by=list(supp,ToothGrowth$dose),FUN=sd)fit<-aov(len~supp*ToothGrowth$dose)summary(fit)detach(ToothGrowth)可視化處理interaction.plot(ToothGrowth$dose,supp,len,type="b",col=c("red","blue"),pch=c(16,18),main="InteractionbetweenDoseandSupplementType")可?gplots包中的plotmans（）來(lái)展?交互效應(yīng)install.packages("gplots")library(gplots)plotmeans(len~interaction(supp,ToothGrowth$dose,sep=""),col=c("red","darkgreen"),connect=list(c(1,3,5),c(2,4,6)),main="InteracitonPlotwith95%CIS",xlab="TreatmentandDoseCombination")?HH包中的interaction2wt（）函數(shù)來(lái)可視化結(jié)果library(HH)interaction2wt(len~supp*ToothGrowth$dose)5、重復(fù)測(cè)量?差分析所謂重復(fù)測(cè)量?差分析，即受試者被測(cè)量不??次。w1b1<-subset(CO2,Treatment=="chilled")w1b1fit<-aov(uptake~conc*Type+Error(Plant/(conc)),w1b1)summary(fit)par(las=2)par(mar=c(10,4,4,2))with(w1b1,interaction.plot(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="Carbondioxideuptakerateumol/m^2sec")數(shù)據(jù)集寬格式（wideformat）：列是變量，?是觀測(cè)值，且???個(gè)受試對(duì)象處理重復(fù)測(cè)量設(shè)計(jì)時(shí)，需要有長(zhǎng)格式（longformat）數(shù)據(jù)才能擬合模型；在長(zhǎng)格式中，因變量每次測(cè)量都要放到它獨(dú)有的?中。reshape包可為?正直將數(shù)據(jù)轉(zhuǎn)換為相應(yīng)的格式。6、多元?差分析（1）單因素多元?差分析library(MASS)head(UScereal)attach(UScereal)y<-cbind(calories,fat,sugars)aggregate(y,by=list(shelf),FUN=mean)cov(y)fit<-manova(y~shelf)summary(fit)summary.aov(fit)（2）評(píng)估假設(shè)檢驗(yàn)單因素多元?差分析有兩個(gè)前提假設(shè)，?個(gè)是多元正態(tài)性，?個(gè)是?差-協(xié)?差同質(zhì)性。前者可?Q-Q圖來(lái)檢驗(yàn)該假設(shè)條件；?差-協(xié)?差矩陣同持性即指各組的協(xié)?差矩陣相同，可?Box'sM檢驗(yàn)來(lái)估計(jì)該假設(shè)。多元正態(tài)分布：若有?個(gè)p*1的多元正態(tài)隨機(jī)向量x，均值為u，存在協(xié)?差矩陣，那么x與u的馬?距離的平?服從?由度為p的卡?分布。Q-Q圖展?卡?頒的分位數(shù)，橫縱坐標(biāo)分別是樣本量與馬?距離平?值。如果點(diǎn)全部落在斜率為1、截距為0的直線上，則表明數(shù)據(jù)服從多元正態(tài)分布。center<-colMeans(y)n<-nrow(y)p<-ncol(y)cov<-cov(y)d<-mahalanobis(y,center,cov)coord<-qqplot(qchisq(ppoints(n),df=p),d,main="Q-QPlotAssessingMultivariateNormality",ylab="MahalanobisD2")abline(a=0,b=1)identify(coord$x,coord$y,labels=s(UScereal))可?mvoutlier包中的ap.plot（）函數(shù)來(lái)檢驗(yàn)多元離群點(diǎn)install.packages("mvoutlier")library(mvoutlier)outliers<-aq.plot(y)outliers（3）穩(wěn)健多元?差分析若多元正態(tài)性或者?差-協(xié)?差均值假設(shè)都不滿?，?擔(dān)?多元離群點(diǎn)，可考慮?穩(wěn)健或?參版本的MANOVA檢驗(yàn)。rrcov包中的Wilks.test（）函數(shù)實(shí)現(xiàn)vegan包中的adonis（）函數(shù)提供了?參數(shù)MANOVA的等同形式library(rrcov)Wilks.test(y,shelf,method="mcd")7、?回歸來(lái)做ANOVA?aov（）函數(shù)擬合模型library(multcomp)levels(cholesterol$trt)fit.aov<-aov(response~trt,data=cholesterol)summary(fit.aov)?回歸lm（）來(lái)解決ANOVA問(wèn)題fit.lm<-lm(response~trt,data=cholesterol)summary(fit.lm)因線性模型要求預(yù)測(cè)變量是數(shù)值型，當(dāng)lm（）函數(shù)碰到因?時(shí)，它會(huì)??系列因??平相對(duì)應(yīng)的數(shù)值型對(duì)照變量為代替因?。若因?有k個(gè)?平，它將會(huì)創(chuàng)建k-1個(gè)對(duì)照變量。contrasts(cholesterol$trt)內(nèi)置對(duì)照組對(duì)照變量創(chuàng)建?法描述contr.helmert第?個(gè)與第?個(gè)?平對(duì)照第三個(gè)?平對(duì)照前兩個(gè)均值第四個(gè)?平對(duì)照前三個(gè)的均值contr.poly基于正交多項(xiàng)式的對(duì)照，?于趨勢(shì)分析和等距?平的有序因?contr.sum對(duì)照變量之和限制為0，也稱(chēng)作偏差找對(duì)，對(duì)各?平的均值與所有?平的