R實(shí)戰(zhàn)學(xué)習(xí)筆記

內(nèi)置數(shù)據(jù)結(jié)構(gòu)：元組：元組由不同元素構(gòu)成，每個元素可以存儲不同類型的數(shù)據(jù)，如字符串、數(shù)值、元組。元組是寫保護(hù)的，元組創(chuàng)建后不能再做任何修改。reload()函數(shù)將以前導(dǎo)入過的模塊再加載一次。重新加載（reload）包括最初導(dǎo)入模塊時應(yīng)用的分析過程和初始化過程。這樣就允許在不退出解釋器的情況下重新加載已更改的Python模塊。Pow(n,2)對某數(shù)求平方pow(4.5-4,2)記（Sqrt()求平方根歐幾里得距離sqrt(pow(4.5-4)+pow(2-1,2))Deepcopy（）深拷貝，能拷貝對象內(nèi)部所有數(shù)據(jù)和引用。若A深拷貝數(shù)據(jù)B，則A變，B不變。Copy（）淺拷貝，只是復(fù)制數(shù)據(jù)。若A淺拷貝數(shù)據(jù)B，則A變，B也變。+表示連接%系統(tǒng)R語言入門R區(qū)分大小寫字母獲取幫助help.start()help(foo)或者？fooexample(foo)data()列出加載包中所有可用示例數(shù)據(jù)集工作空間工作空間是R用來讀取文件和保存結(jié)果的默認(rèn)目錄。getwd（）setwd(“D:/rworkspace”)D:/rworkspace必須存或者用dir.create()創(chuàng)建新目錄load（）從上一次會話結(jié)束地方開始輸入source（”myscript.R”）將執(zhí)行包含在文件myscript.R中的R語句集合輸出sink(“myoutput”,append=TRUE,split=TRUE)默認(rèn)情況下，有一樣名字文件，會被覆蓋原文件。append將文本追加到文件后，而不是覆蓋，split將文件發(fā)送至屏幕并且輸出文件。不加參數(shù)的sink（）僅向屏幕輸出結(jié)果。圖形輸出Pdf(“mygraphs.pdf”)輸出結(jié)果作為輸入包包是函數(shù)、數(shù)據(jù)、預(yù)編譯代碼以一種定義完善的格式組成的集合。存包的目錄稱為庫函數(shù).libpath()顯示庫所在位置library()顯示庫位置以及有哪些包第一次安裝一個包用:install.packages(“gclus”)包（引號是必須的）gclus提供了創(chuàng)建增強(qiáng)型散點(diǎn)圖的函數(shù)。包只需安裝一次，但作者可能更新它update.packages(“包名”)查看已經(jīng)安裝的包installed.packages()列出安裝包的版本和依賴關(guān)系包的載入：包安裝時從cran鏡像站點(diǎn)下載它加入庫的過程需用到時還得library（gclus）包帶有示例，help（package=”包名”）可以輸出包的簡短描述以及包中函數(shù)和數(shù)據(jù)集名稱的列表。選中全部代碼再按run否則編譯器只會運(yùn)行最后一行代碼數(shù)據(jù)集對象是可以賦值給變量的任何事物。x<-mtcars[order(mpg),]不加,系統(tǒng)分不清向量矩陣數(shù)組數(shù)據(jù)框data.frame()patientid$age用來選取一個給定數(shù)據(jù)框中的某個特定變量attach(mtcars)//可以將數(shù)據(jù)框添加到R的搜索路徑，plot(disp,mpg)//此時若有其他名稱相同對象，比如也叫disp則原始對象取得優(yōu)先權(quán)。若長度不匹配則報(bào)錯，要注意detach（mtcars）否則得plot(mtcars$disp,mtcars$mpg)with（{只寫半邊括號，r中命令就能按enter后自動換行命令<<-按上下鍵可以使用歷史命令列表（向量、數(shù)組、數(shù)據(jù)框、列表）因子：因子(factors)提供了一種處理分類數(shù)據(jù)的更簡介的方式。levels=c(“poor”,””,”improved”,”excellent”)可以通過level將類別變量排序列聯(lián)表table(patientdata$diabete,patientdata$status)str(patientdata)顯示數(shù)據(jù)框patientdata的結(jié)構(gòu)，如哪些列，多少個變量及觀測，具體值是。Summary（patient）顯示統(tǒng)計(jì)概要#注釋無多行注釋圖形初階如何創(chuàng)建多個圖形并隨時查看每一個呢方法一dev.new()創(chuàng)建圖語句dev.new()創(chuàng)建圖語句dev.off()將輸出返回到終端修改圖形參數(shù)>dose<-c(20,30,40,45,60)>drugA<-c(16,20,27,40,60)>drugB<-c(15,18,25,31,40)>plot(dose,drugA,type="b")復(fù)制當(dāng)前圖形參數(shù)opar<-par(no.readonly=TRUE)par(lty=2,pch=17)#設(shè)置圖形參數(shù)，此時應(yīng)用于所有之后創(chuàng)建圖形plot(dose,drugA,type="b")#查看修改參數(shù)后圖片效果>par(opar)>plot(dose,drugA,type="b")#圖片又變回原來參數(shù)了>plot(dose,drugA,type="b",lty=2,pch=17)#僅應(yīng)用于該圖像繪制點(diǎn)的符號R中查看已經(jīng)安裝的包路徑和包名library()pchAxis坐標(biāo)軸font字體類型：1常規(guī)2粗體3斜體4粗斜體圖形前景色及背景色圖形尺寸：代碼應(yīng)該提前，不要等到圖已經(jīng)做出來了再寫該代碼Par(pin=c(4,3),m,mar=c(1,.5,1,.2))pin長寬Mar指邊界大小par(mar=(5,4,4,8)+0.1)下開始，逆時針轉(zhuǎn)注意：某些高級繪圖函數(shù)已經(jīng)包含了默認(rèn)的標(biāo)題和標(biāo)簽。你可以通過在plot（）語句或者單獨(dú)的par（）中添加ann=FALSE來移除它們announcement通告坐標(biāo)軸yaxt="n"xaxt=”n”將分別禁用x軸、y軸(會留下框架，只是禁用了刻度)參考線Abline圖例Legend文本標(biāo)注Textmtexttext(wt,mpg,s(mtcars),cex=.6,pos=4,col="red")直接點(diǎn)的位置，再加上標(biāo)注向量text(3,3,”exzample”)單獨(dú)點(diǎn)加標(biāo)注圖形組合Par(mfrow=c(2,2))2行2列按行填充4個圖layout(matrix(c(1,1,2,3),2,2,byrow=TRUE),widths=c(3,1),heights=c(1,2))c()為矩陣元素2,2為2行2列l(wèi)ayout(matrix(c(1,2,3,0,2,3,0,0,3),nr=3))matrix有9個元素，具有這樣的形式：

[,1][,2][,3][1,]

1

0

0[2,]

2

2

0[3,]

3

3

3把這個矩陣傳入layout函數(shù)，我們就能得到這樣的outputdevicepar(fig=c(0,0.8,0,0.8))>plot(wt,mpg,xlab="mile",ylab="weight")par(fig=c(0,0.8,0.55,1),new=TRUE)boxplot(wt,horizontal=TRUE,axes=FALSE)opar<-par(no.readonly=TRUE)par(fig=c(0,0.8,0,0.8))plot(wt,mpg,xlab="mile",ylab="weight")par(fig=c(0,0.8,0.55,1),new=TRUE)boxplot(wt,horizontal=TRUE,axes=FALSE)六、基本圖形：條形圖barplot()一個數(shù)值堆砌條形圖分組條形圖均值條形圖可以使用數(shù)據(jù)整合函數(shù)并將結(jié)果傳遞給barplot()條形圖的微調(diào)棘狀圖棘狀圖對堆砌條形圖進(jìn)行了重縮放，這樣每個條形的高度均為1，每一段高度表示比例。library(vcd)attach(Arthritis)counts<-table(Treatment,Improved)spine(counts,main="SpinogramExample")#spine屬于vcd包中detach(Arthritis)餅圖opar<-par(no.readonly=TRUE)par(mfrow=c(2,2))x<-c(10,12.9,4,16,8)lbls<-c("US","France","Germany","Australia","UK")pie(x,lbls,main="簡單餅圖")lbls<-paste(lbls,"",round(x/sum(x)*100,2),"%",sep="")#paste()實(shí)現(xiàn)字符串連接pie(x,lbls,main="使用百分比標(biāo)簽的餅圖",col=rainbow(length(lbls)))install.packages("plotrix")library(plotrix)pie3D(x,labels=lbls,explode=0.1,main="3D餅圖")mytable<-table(state.region)lbls<-paste(names(mytable),"\n",mytable,sep="")pie(mytable,labels=lbls,main="列表數(shù)據(jù)餅圖\n(一個簡單的)")扇形圖比較各個部分大小par(mfrow=c(1,1))

lbls

<-

c("US",

"France",

"Germany",

"Australia",

"UK")#lblslabellist

fan.plot(x,

labels=lbls,

main="扇形圖")

直方圖histogrambreaks=12,表示分組數(shù)hist（x）hist(mtcars$mpg,breaks=12,col="red",xlab="milespergallon",main="coloredhistogramwith12bins")分成12組#根據(jù)概率密度畫圖#增加軸須圖hist(mpg,freq=FALSE,#fre=FALSE表示根據(jù)概率密度而不是頻數(shù)繪制圖形breaks=12,##

breaks用于控制組的數(shù)量col="red",main="直方圖、軸須圖、概率密度曲線",xlab="MilesPerGallon",ylab="概率密度")rug(jitter(mpg))增加軸須圖#rug(jitter(mpg,amount=0.01))如果數(shù)據(jù)中有許多結(jié)（數(shù)據(jù)中相同的值），可以使用該代碼將軸須圖的數(shù)據(jù)打散lines(density(mpg),col="blue",lwd=2)#概率密度曲線核密度曲線軸須圖：是實(shí)際數(shù)據(jù)值的一種一維呈現(xiàn)方式x<-mpgh<-hist(x,breaks=12,col="red",main="添加正態(tài)密度曲線和外框的直方圖",xlab="MilesPerGallon",ylab="頻數(shù)")xfit<-seq(min(x),max(x),length=40)#生成一個序列，最小最大值為···長度為40個yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))#生成一個yfit<-yfit*diff(h$mids[1:2])*length(x)不理解lines(xfit,yfit,col="blue",lwd=2)box()#外框detach(mtcars)par(opar)核密度圖核密度估計(jì)是用于估計(jì)隨機(jī)變量概率密度函數(shù)的一種非參數(shù)方法。Plot(density(x))可比較核密度圖箱線圖attach(mtcars)boxplot(mpg,main="Boxplot",ylab="milespergallon")并列箱線圖跨組比較boxplot(formula,data=data.frame)formula為一個公式，Y~A為類別型變量A每個值并列生成數(shù)值型變量y的箱線圖Y~A*B為類別型變量A和B所有水平的兩兩組合生成數(shù)值型變量y的箱線圖varWidth=TRUE箱線圖寬度及樣本大小的平方根成正比。boxplot(mpg~cyl,data=mtcars,main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")boxplot(mpg~cyl,data=mtcars,notch=TRUE,varwidth=TRUE,main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")notch=TRUE得到含凹槽的箱線圖，若兩個箱的凹槽互不重疊則表明他們的中位數(shù)顯著差異。兩個交叉因子的箱線圖cyl.f氣缸468am.f,自動及標(biāo)準(zhǔn)型號cyl.f<-factor(cyl,levels=c(4,6,8),labels=c("4clinder","6cylinder","8cylinder"))>am.f<-factor(am,levels=c(0,1),labels=c("auto","standard"))>boxplot(mpg~cyl.f*am.f,data=mtcars,varwidth=TRUE,col=c("gold","darkgreen"),main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")am.f*cyl.f這個順序影響排序小提起圖是箱線圖和核密度圖的合體‘散點(diǎn)圖Plot（）點(diǎn)圖dotchart()attach(mtcars)x<-mtcars[order(mpg),]x$cyl<-factor(x$cyl)hex$color[x$cyl==4]<-"red"x$color[x$cyl==6]<-"blue"x$color[x$cyl==8]<-"darkgreen"dotchart(x$mpg,labels=s(x),cex=.7,groups=x$cyl,gcolor="black",color=x$color,pch=19,main="gasmileageforcarmodels\ngroupedbycylinder",xlab="milespergallon")groups=x$cyl,分組x$color字符型向量color被添加到數(shù)據(jù)框中創(chuàng)建分組因子factor(x=character(),levels,labels=levels,exclude=NA,ordered=is.ordered(x),nmax=NA)七、描述性統(tǒng)計(jì)分析總體描述library(psych)attach(mtcars)vars<-c("mpg","hp","wt")describe(mtcars[vars])detach(mtcars)結(jié)果分組描述性統(tǒng)計(jì)量library(reshape)dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))dfm<-melt(mtcars,measure.vars=c("mpg","hp","wt"),id.vars=c("am","cyl"))cast(dfm,am+cyl+variable~.,dstats)#黃色為類別型分組變量，右邊為特征統(tǒng)計(jì)量結(jié)果列聯(lián)表考察交叉類型的頻數(shù)、比例、以及交叉類型之間的相關(guān)性一維>mytable<-with(Arthritis,table(Improved))>mytableImprovedNoneSomeMarked421428二維install.packages("gtools")install.packages("gdata")install.packages("gmodels")library(gmodels)with(Arthritis,CrossTable(Treatment,Improved))Help(crosstable)可以計(jì)算小數(shù)位數(shù)、卡方等獨(dú)立性檢驗(yàn)、計(jì)算期望皮爾遜殘差等。三維其他prop.table()計(jì)算比例Margin.table()計(jì)算頻數(shù)獨(dú)立性、相關(guān)性考察數(shù)據(jù)分析一文件八、回歸0LS（ordinaryleastsquare）普通最小二乘法的回歸簡單線性回歸殘差標(biāo)準(zhǔn)誤可以理解為:模型用身高預(yù)測體重的平均誤差。R平方：=0.991表明模型可以解釋99.1%的方差，它也是實(shí)際及預(yù)測值之間的相關(guān)系數(shù)。F統(tǒng)計(jì)量：等同于身高回歸系數(shù)的t檢驗(yàn)分析所以可以用一個彎曲曲線來提高預(yù)測的精度。多項(xiàng)式回歸P<0.001水平下，回歸系數(shù)均非常顯著。模型方差解釋率達(dá)到99.9%。二次項(xiàng)顯著性（t=13.891,p<0.001）表明二次項(xiàng)提高了模型的擬合度。（這樣就可以拒絕原假設(shè)二次項(xiàng)系數(shù)為0的假設(shè)）擬合曲線Weight=261.88-7.35heght+0.083height^2多元線性回歸預(yù)測變量不止一個。例如考察人口、文盲率、收入、霜凍綜合對謀殺率的影響。案例：檢測變量關(guān)系，相關(guān)系數(shù)矩陣>states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])>cor(states)MurderPopulationIlliteracyIncomeFrostMurder1.00000000.34364280.7029752-0.2300776-0.5388834Population0.34364281.00000000.10762240.2082276-0.3321525Illiteracy0.70297520.10762241.0000000-0.4370752-0.6719470Income-0.23007760.2082276-0.43707521.00000000.2262822Frost-0.5388834-0.3321525-0.67194700.22628221.0000000做散點(diǎn)圖矩陣>library(car)scatterplotMatrix(states,spread=FALSE,lty.smooth=2,main="scatterplotmatrix")從圖中可以看出謀殺率是雙峰曲線，謀殺率隨著人口文盲率增加而增加，隨著收入水平和結(jié)霜天數(shù)增加而下降。越冷地方文盲率越低，收入水平越高。多元線性回歸參數(shù)查看DFdegreeoffreedom自由度有交互項(xiàng)的多元線性回歸許多有趣的研究都會涉及交互的預(yù)測變量。見書R實(shí)戰(zhàn)170頁用圖形展示上圖藍(lán)色文字含義回歸診斷對模型的信賴程度依賴于它在多大程度上滿足OLS模型的統(tǒng)計(jì)假設(shè)。OLS回歸的統(tǒng)計(jì)假設(shè)正態(tài)性:QQ圖獨(dú)立性：因變量是否獨(dú)立。例如Weight~height。一位女性體重會影響另一位女性體重？如果來自同一個家庭有可能線性同方差性模型在多大程度上滿足統(tǒng)計(jì)假設(shè)的任何信息。標(biāo)準(zhǔn)方法#8.3.1回歸診斷標(biāo)準(zhǔn)方法fit<-lm(weight~height,data=women)par(mfrow=c(2,2))plot(fit)上圖沒有顯示因變量是否獨(dú)立的圖改進(jìn)的方法進(jìn)行回歸診斷Car包提供了大量函數(shù)，增強(qiáng)了對擬合回歸模型評價的能力。正態(tài)性方法一：檢驗(yàn)正態(tài)性par(mfrow=c(1,1))library(car)states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)qqPlot(fit,labels=s(states),id.method="identify",simulate=TRUE,main="Q-QPlot")#id.method="identify'能夠交互式繪圖，鼠標(biāo)單擊圖形內(nèi)的點(diǎn)，將會標(biāo)注函數(shù)中l(wèi)abel選項(xiàng)的設(shè)定值。#simulate=TRUE95%的置信區(qū)間將會用參數(shù)自助法關(guān)注下離群點(diǎn)方法二：殘差圖residplot<-function(fit,nbreaks=10){z<-rstudent(fit)hist(z,breaks=nbreaks,freq=FALSE,xlab="StudentizedResidual",main="DistributionofErrors")rug(jitter(z),col="brown")curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=2)lines(density(z)$x,density(z)$y,不好理解col="red",lwd=2,lty=2)legend("topright",legend=c("NormalCurve","KernelDensityCurve"),lty=1:2,col=c("blue","red"),cex=.7)residplot(fit)方法三：標(biāo)準(zhǔn)方法中的QQ圖，看是否在一條45度直線上。誤差的獨(dú)立性因變量值（或殘差）是否相互獨(dú)立，最好依據(jù)收集數(shù)據(jù)方式的先驗(yàn)知識。比如：時間序列呈現(xiàn)自相關(guān)性，相隔時間越近的觀測相關(guān)性大于相隔越遠(yuǎn)的觀測。Car包中Durbin-Watson函數(shù)作Durbin-Watson檢驗(yàn)

線性通過成分殘差圖即偏殘差圖，可以看看因變量及自變量之間是否線性相關(guān)。library(car)crPlots(fit)綠色線為平滑擬合曲線。同方差性Car包中提供了ncvTest()和spreadLevelPlot（）函數(shù)判斷方差是否恒定。ncvTest()函數(shù)生成一個計(jì)分檢驗(yàn)，零假設(shè)為誤差方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化。若檢驗(yàn)顯著，則說明存在異方差性，即誤差方差不恒定。spreadLevelPlot（）創(chuàng)建一個添加了最佳擬合曲線的散點(diǎn)圖。展示標(biāo)準(zhǔn)化殘差絕對值及擬合值的關(guān)系。spreadLevelPlot(fit)線性模型假設(shè)的綜合驗(yàn)證library(gvlma)gvmodel<-gvlma(fit)summary(gvmodel)多重共線性回歸系數(shù)測量的是當(dāng)其他預(yù)測變量不變時，某個預(yù)測變量對響應(yīng)變量的影響。比如：假定年齡不變，測量握力及年齡的關(guān)系？。這種問題稱為多重共線性。在回歸模型牽涉到多個自變量的時候，自變量之間可能會相互關(guān)聯(lián)，即他們之間存在有多重共線性。它會導(dǎo)致模型參數(shù)的置信區(qū)間過大，使得單個系數(shù)解釋起來很困難。方差膨脹因子：VIF（varianceinflationfactor）VIF的平方根表示變量回歸參數(shù)的置信區(qū)間能膨脹為及模型無關(guān)的預(yù)測變量的程度。異常觀測值全面回歸分析需要包含對異常值的分析。離群點(diǎn)模型預(yù)測效果不佳的觀測點(diǎn)，殘差比較大的點(diǎn)。判斷：1：Q—Q圖中落在置信區(qū)間以外的點(diǎn)被認(rèn)為是離群點(diǎn)2：標(biāo)準(zhǔn)化殘差值大于2或者小于-2可能是離群點(diǎn)Car包中outlierTest()函數(shù)可以求得最大標(biāo)準(zhǔn)化殘差絕對值Bonferroni調(diào)整后的P值。P=0.048該函數(shù)只是根據(jù)單個最大殘差值得顯著性判斷是否有離群點(diǎn)。若不顯著則沒有離群點(diǎn)，若顯著則需要刪除離群點(diǎn)，再檢驗(yàn)是否有其他離群點(diǎn)存在高杠桿值點(diǎn)？及響應(yīng)變量值沒有關(guān)系，由許多異常的預(yù)測變量值組合起來的。通過帽子統(tǒng)計(jì)量判斷。hat.plot<-function(fit){p<-length(coefficients(fit))包括截距項(xiàng)若四個預(yù)測變量則為5n<-length(fitted(fit))觀測數(shù)量，比如有50個樣本值plot(hatvalues(fit),main="IndexPlotofHatValues")abline(h=c(2,3)*p/n,col="red",lty=2)#水平線為帽子均值2倍和3倍的位置。identify(1:n,hatvalues(fit),names(hatvalues(fit)))hat.plot(fit)p/n即為帽子均值。lm(Murder~Population+Illiteracy+Income+Frost,data=states)謀殺率大小由4個預(yù)測變量及一個截距（常數(shù)項(xiàng)）分別貢獻(xiàn)。強(qiáng)影響點(diǎn)對模型參數(shù)估計(jì)值影響有些比例失衡的點(diǎn)。比如若移除模型的一個觀測點(diǎn)時，模型會發(fā)生巨大的變化，則需要檢測數(shù)據(jù)中是否存在強(qiáng)影響點(diǎn)。方法：cook距離（D統(tǒng)計(jì)量）、CookD值大于4/(n-k-1)則表明她是強(qiáng)影響點(diǎn)。N為樣本量大小，k為預(yù)測變量的數(shù)目。cutoff

<-

4/(nrow(states)-length(fit$coefficients)-2)

plot(fit,

which=4,

cook.levels=cutoff)

abline(h=cutoff,

lty=2,

col="red")方法2：變量添加圖：上圖可以鑒別強(qiáng)影響點(diǎn)，但是并不提供這些點(diǎn)如何影響模型的參數(shù)的信息。library(car)avPlots(fit,ask=FALSE,onepage=TRUE,id.method="identify")主要看紅線2邊是否只有一邊有點(diǎn)，另外一邊沒有點(diǎn)。將離群點(diǎn)、杠桿值和強(qiáng)影響點(diǎn)整合在一張圖畫中：influencePlot(fit,id.method="identify",main="InfluencePlot",sub="CirclesizeisproportionaltoCook'sdistance")id.method="identify"為交互作圖，點(diǎn)擊感興趣的點(diǎn)，即顯示信息。對離群點(diǎn)的改進(jìn)措施如何處理違背回歸假設(shè)的問題。刪除觀測點(diǎn)如果推斷離群點(diǎn)是數(shù)據(jù)記錄錯誤等，可以選擇將最大離群點(diǎn)、強(qiáng)影響點(diǎn)刪除。刪除后重新擬合，若仍然有離群點(diǎn)及強(qiáng)影響點(diǎn)，則重復(fù)上面步驟，直至獲得滿意結(jié)果。若不了解產(chǎn)生離群點(diǎn)原因，則需謹(jǐn)慎對待，不要盲目刪除。變量變換當(dāng)模型不符合正態(tài)性、線性、或者同方差性假設(shè)時，可以將一個或者多個變量變換。變換多用Yλ代替Yλ常用取值-2

-1

-0.5

0

0.5

1

2若Y是比例數(shù)，通常使用logit變換[ln(Y/1-Y)]當(dāng)違反了正態(tài)性假設(shè)常常對響應(yīng)變量Y嘗試某種變換。Box-cox正態(tài)變換當(dāng)違反了線性假設(shè)時，對預(yù)測變量X進(jìn)行變換常常會比較有用。當(dāng)違反方差穩(wěn)定性，可以變換響應(yīng)變量。具體見上文spreadLevelPlot()謹(jǐn)慎對待變量變換：如果變換了變量，你的解釋必須機(jī)遇變換后的變量，而不是初始變量。如果變換得沒有意義，則應(yīng)該避免這樣做。比如：收入的對數(shù)變換好理解。但是自殺意念的頻率域抑郁癥程度的立方根之間的關(guān)系在現(xiàn)實(shí)中很難解釋出有意義的含義。添加或者刪除變量是處理多重共線性的非常重要的方法。僅僅做預(yù)測則多重共線性并不構(gòu)成問題，但是如果要解釋每個預(yù)測變量，則必須解決該問題。常用方法：刪除某個存在多重共線性的變量（vif>2嶺回歸—多元回歸的變體使用其他回歸模型如果存在離群點(diǎn)或者強(qiáng)影響點(diǎn)可以使用穩(wěn)健回歸模型替代OLS回歸模型。如果違背正態(tài)性假設(shè)可以用非參數(shù)回歸模型。如果存在顯著非線性，嘗試用非線性回歸模型。如果違背了誤差獨(dú)立性準(zhǔn)則，嘗試用專門研究誤差結(jié)構(gòu)的模型，如：時間序列模型，或者多層次回歸模型，或者廣義線性模型。選擇“最佳”回歸模型預(yù)測精度及模型簡潔度的調(diào)和問題。最佳是一種權(quán)衡，沒有絕對最佳模型。模型比較比較兩個嵌套模型的擬合優(yōu)度。嵌套式指該模型的一些項(xiàng)，完全包含在另一個模型中。在States多元回歸模型中可以發(fā)現(xiàn)Income和frost回歸系數(shù)不顯著，檢驗(yàn)不含這兩個變量的模型預(yù)測效果是否更好。方法一：anova()函數(shù)方法二：AIC赤池信息準(zhǔn)則ANOVA需要嵌套模型，而AIC方法不需要變量選擇從大量變量中選擇最終的預(yù)測變量逐步回歸模型會一次添加或者刪除一個變量，直到達(dá)到某個判停準(zhǔn)則為止。向前逐步回歸：每次添加一個變量向后逐步回歸：每次刪除一個變量向前向后逐步回歸（簡稱逐步回歸）：變量每次進(jìn)入一個，對貢獻(xiàn)不大變量再刪除。逐步回歸也許找到個好模型，但是不能保證模型就是最佳模型，因?yàn)椴皇敲恳粋€模型都被評價了，故有了全子集回歸法。Intercept截距項(xiàng)疑問：AIC怎么計(jì)算的？全子集回歸所有可能的模型都會被檢驗(yàn)。全子集回歸可用leaps包中regsubsets（）函數(shù)實(shí)現(xiàn)，通過R平方，調(diào)整R平方或者mallowscp統(tǒng)計(jì)量等準(zhǔn)則來選擇最佳模型。R平方含義是預(yù)測變量解釋響應(yīng)變量的程度。一般R平方會隨著預(yù)測變量的增加而增加，但是當(dāng)及樣本量相比，預(yù)測變量數(shù)目很大時，容易導(dǎo)致過擬合。R平方很可能會丟失數(shù)據(jù)的偶然變異信息，而調(diào)整的R平方提供了更為真實(shí)的R平方估計(jì)。install.packages("leaps")

library(leaps)

leaps

<-

regsubsets(Murder

~

Population

+

Illiteracy

+

Income

+

Frost,

data=states,

nbest=4)

plot(leaps,

scale="adjr2")mallowscp統(tǒng)計(jì)量也是作為逐步回歸的判停準(zhǔn)則。好的回歸模型它的cp統(tǒng)計(jì)量非常接近模型的參數(shù)數(shù)目（包含截距項(xiàng)）library(car)

subsets(leaps,

statistic="cp",

main="Cp

Plot

for

All

Subsets

Regression")

abline(1,

1,

lty=2,

col="red")注意：該圖也是交互作圖，有時運(yùn)行沒有停止運(yùn)行，點(diǎn)擊下圖看是否為大部分情況全子集回歸優(yōu)于逐步回歸，但是數(shù)據(jù)量大時全子集回歸會很慢。擬合效果佳但是沒有意義的模型沒有意義。主題背景知識的理解才能指導(dǎo)獲得理想的模型。深層分析模型泛化能力和相對重要性交叉驗(yàn)證一定比例數(shù)據(jù)作訓(xùn)練集，一定比例數(shù)據(jù)作檢驗(yàn)?zāi)Ｐ蜆颖?。K重交叉驗(yàn)證：樣本被分為K個子樣本，輪流將K-1個子樣本組合為訓(xùn)練集，另外一個子樣本作為保留集合。得到K個預(yù)測方程，記錄K個保留樣本的預(yù)測表現(xiàn)結(jié)果，然后求其平均值。#

代碼清單8-15

R平方的k重交叉驗(yàn)證

shrinkage

<-

function(fit,

k=10){

require(bootstrap)

theta.fit

<-

function(x,y){lsfit(x,y)}

theta.predict

<-

function(fit,x){cbind(1,x)%*%fit$coef}

x

<-

fit$model[,2:ncol(fit$model)]

y

<-

fit$model[,1]

results

<-

crossval(x,

y,

theta.fit,

theta.predict,

ngroup=k)

r2

<-

cor(y,

fit$fitted.values)^2

r2cv

<-

cor(y,

results$cv.fit)^2

cat("Original

R-square

=",

r2,

"\n")

cat(k,

"Fold

Cross-Validated

R-square

=",

r2cv,

"\n")

cat("Change

=",

r2-r2cv,

"\n")

}#K重交叉驗(yàn)證函數(shù)

fit1

<-

lm(Murder

~

Population

+

Income

+

Illiteracy

+

Frost,

data=states)

shrinkage(fit1)

fit2

<-

lm(Murder

~

Population

+

Illiteracy

,

data=states)

shrinkage(fit2)相對重要性哪些變量對預(yù)測最為重要。根據(jù)相對重要性對預(yù)測變量排序。方法一：預(yù)測變量若不相關(guān)，則根據(jù)預(yù)測變量及響應(yīng)變量的相關(guān)系數(shù)即可確認(rèn)重要性排序。但大部分情況下預(yù)測變量之間是有一定相關(guān)性。方法二標(biāo)準(zhǔn)化的回歸系數(shù)，表示當(dāng)其他預(yù)測變量不變時，該預(yù)測變量一個標(biāo)準(zhǔn)差的變化可引起的響應(yīng)變量的預(yù)期變化。方法三：相對權(quán)重，測量對R平方的貢獻(xiàn)附錄：#

代碼清單8-16

relweights()函數(shù)，計(jì)算預(yù)測變量的相對權(quán)重

relweights

<-

function(fit,

...){

R

<-

cor(fit$model)

nvar

<-

ncol(R)

rxx

<-

R[2:nvar,

2:nvar]

rxy

<-

R[2:nvar,

1]

svd

<-

eigen(rxx)

evec

<-

svd$vectors

ev

<-

svd$values

delta

<-

diag(sqrt(ev))

lambda

<-

evec

%*%

delta

%*%

t(evec)

lambdasq

<-

lambda^2

beta

<-

solve(lambda)

%*%

rxy

rsquare

<-

colSums(beta^2)

rawwgt

<-

lambdasq

%*%

beta^2

import

<-

(rawwgt/rsquare)

*

100

lbls

<-

names(fit$model[2:nvar])

rownames(import)

<-

lbls

colnames(import)

<-

"Weights"

barplot(t(import),

names.arg=lbls,

ylab="%

of

R-Square",

xlab="Predictor

Variables",

main="Relative

Importance

of

Predictor

Variables",

sub=paste("R-Square=",

round(rsquare,

digits=3)),

...)

return(import)

}九、方差分析基本術(shù)語ANOVA：(AnalysisofVariance方差分析:關(guān)注組別差異的分析。組間因子：每個患者僅被分配到一個組別中。組內(nèi)因子：如植物類型和二氧化碳濃度。每個植物只能劃分到一個植物類型A或者B中。但是所A和B中所有植物均可以在相同濃度下測試二氧化碳吸收量，則二氧化碳濃度為組內(nèi)因子。均衡設(shè)計(jì)：各個組觀測數(shù)目相同組內(nèi)因子：時間是兩水平（五周、六周）的組內(nèi)因子，單因素組內(nèi)方差分析：每個患者在所有水平下都進(jìn)行了測量。重復(fù)測量方差分析：每個受試者不止一次被測量。多元方差分析：因變量不止一個多元協(xié)方差分析：協(xié)變量存在。協(xié)變量：白氏抑郁癥量表（BDI）記錄了他們的抑郁水平，那么你可以評測療法類型的影響前，對任何抑郁水平的組間差異進(jìn)行統(tǒng)計(jì)性調(diào)整。協(xié)方差分析：BDI稱為協(xié)變量，該設(shè)計(jì)為協(xié)方差分析。ANOVA（）模型擬合ANOVA和回歸方法都是獨(dú)立發(fā)展而來，但是從函數(shù)形式上看，他們都是廣義線性模型的特例。lm()函數(shù)也能分析ANOVA模。但是一般用aov()函數(shù)展示結(jié)果。aov(formula，data=dataframe)表達(dá)式中各項(xiàng)的順序因子不止一個，并且是非平衡設(shè)計(jì)；存在協(xié)變量。出現(xiàn)任意一種情況，等式右邊的變量都及其他變量相關(guān)，此時，無法清晰的劃分他們對因變量的影響。表達(dá)式的順序很重要#

y

~

A

#

y

~

x

+

A

#

y

~

A

*

B

#

y

~

x1

+

x2

+

A*B

#

y

~

B

+

A

#

y

~

A

+

Error(Subject/A)

#

y

~

B

*

W

+

Error(Subject/W)對雙因素方差分析，若觀測數(shù)不同，則Y~A*B及Y~B*A的結(jié)果不同。1序貫型Y~A+B+A：BR中的ANOVA表的結(jié)果將評價：A對Y的影響控制A時，B對Y的影響。控制A和B的主效應(yīng)時，A及B的交互效應(yīng)。2分層型：A根據(jù)B調(diào)整，B依據(jù)A調(diào)整，A:B交互項(xiàng)同時根據(jù)A和B調(diào)整3邊界型：A根據(jù)B和A:B做調(diào)整，A:B交互項(xiàng)根據(jù)A和B調(diào)整樣本大小越不平衡，效應(yīng)項(xiàng)的順序?qū)Y(jié)果的影響越大。一般來說越基礎(chǔ)的效應(yīng)越需要放在表達(dá)式前面。具體而言，首先是協(xié)變量，然后是主效應(yīng)，接著是雙因素的交互項(xiàng)，再接著是三因素的交互項(xiàng)，以此類推。對主效應(yīng)，越基礎(chǔ)性的變量越應(yīng)該放在表達(dá)式前面，因此性別放在處理方式之前。單因素方差分析比較分類因子定義的兩個或者多個組別中因變量的均值。例子：50個患者接受5種治療方案，求哪種藥物療法降低膽固醇（響應(yīng)變量）最多繪制帶有置信區(qū)間的均值圖形。繪制帶有95%置信區(qū)間的各療法均值。>library(gplots)>plotmeans(response~trt,xlab="Treatment",ylab="Response",main="MeanPlot\nwith95%CI")>detach(cholesterol)95%置信區(qū)間的各療法均值藥物E降低膽固醇最多，而1time降低最少多重比較F檢驗(yàn)表明5種藥物療效不同，但是并沒有說明哪種療效及其他療效不同。TukeyHSD(fit)

par(las=2)

par(mar=c(5,8,4,2))

plot(TukeyHSD(fit))疑問：用什么函數(shù)得出組間差異置信區(qū)間的，Ｆ檢驗(yàn)？方法二：#有相同字母的組說明均值差異不顯著library(multcomp)par(mar=c(5,4,6,2))tuk<-glht(fit,linfct=mcp(trt="Tukey"))plot(cld(tuk,level=.05),col="lightgrey")評估檢驗(yàn)的假設(shè)條件單因素方差分析結(jié)果準(zhǔn)確性，依賴于可做單因素方差分析的假設(shè)。假設(shè)有假設(shè)因變量服從正態(tài)分布各組方差相等具體檢驗(yàn)方法假設(shè)因變量服從正態(tài)分布用QQ圖檢驗(yàn)正態(tài)性各組方差相等用方差齊次性檢驗(yàn)函數(shù)（對離群點(diǎn)敏感）用方差齊次性檢驗(yàn)函數(shù)有Bartlett檢驗(yàn)、Fligner-Killeen檢驗(yàn)、Brown-Forsythe檢驗(yàn)（HH包中hov()函數(shù)）。選一個就行，結(jié)果相同。要看看是否有離群點(diǎn)使得各組間方差沒有顯著不同故真的是各組方差數(shù)據(jù)沒有顯著差異。結(jié)論：該數(shù)據(jù)可以用ANOVA模型擬合得非常好。單因素協(xié)方差分析包含了一個或者多個定量的協(xié)變量。舉個例子就是：降雨量（t）=K*溫度（t）+e其中，t是自變量時間，降雨量（t）是因變量，而溫度（t）則是協(xié)變量K為一個常數(shù)。例子：自變量（dose）懷孕小組分為四組，每組接受（0,5,50或者500）的藥物處理。產(chǎn)下幼仔重量均值為因變量，懷孕時間為協(xié)變量。例子：自變量（dose）懷孕小組分為四組，每組接受（0,5,50或者500）的藥物處理。產(chǎn)下幼仔重量均值為因變量，懷孕時間為協(xié)變量。通過F檢驗(yàn)，檢驗(yàn)兩組的方差是否有顯著性差異。S^2=∑(X-X平均)^2/(n-1)兩組數(shù)據(jù)就能得到兩個S^2值，S大^2和S小^2F=S大^2/S小^2由表中f大和f?。╢為自由度n-1),查得F表，然后計(jì)算的F值及查表得到的F表值比較，如果F<F表表明兩組數(shù)據(jù)沒有顯著差異；F≥F表表明兩組數(shù)據(jù)存在顯著差異置信度95%時F值(單邊)f大f小2345678910∞2345678910∞19.09.556.945.795.144.744.464.264.103.0019.169.286.595.414.764.354.073.863.713.6019.259.126.395.194.534.123.843.633.482.3719.309.016.265.054.393.973.693.483.333.2119.338.946.164.954.283.873.583.373.222.1019.368.886.094.884.213.793.503.293.142.0119.378.846.044.824.513.733.443.233.071.9419.388.816.004.784.103.683.393.183.021.8819.398.785.964.744.063.633.343.132.971.8319.58.535.634.363.673.232.932.712.541.00獲取調(diào)整的各組均值，即去除協(xié)變量效應(yīng)后各組均值。F檢驗(yàn)表明了不同劑量下各幼鼠體重均值不同，但是沒有告訴我們哪種處理方式及其他方式不同。（此處是下多大劑量及其他組不同）T值為正值（2.581）說明未用藥比用藥體重高。其他對照調(diào)節(jié)rbind()詳見help(glht)GeneralLinearHypotheses一般性線性假設(shè)平估檢驗(yàn)的假設(shè)條件假設(shè)正態(tài)性同方差性回歸斜率相同：本例中假定四個處理組通過懷孕時間（協(xié)變量）來預(yù)測出生體重（因變量）的回歸斜率相同。自變量一般為類別型變量。（協(xié)變量定義時就有這一點(diǎn)，協(xié)變量及因變量必須呈現(xiàn)線性關(guān)系）ANCOVA模型包含懷孕時間*劑量的交互項(xiàng)時如：（Weight=gesttime+dose+gesttime:dose），可對回歸斜率的同質(zhì)性進(jìn)行檢驗(yàn)。交互效應(yīng)若顯著，則意味著懷孕時間和幼崽出生體重間的關(guān)系依賴于藥物劑量的水平。當(dāng)然不顯著才行···交互不顯著，則模型約為：Weight=gesttime+dose此時回歸斜率就基本相同，只是截距不同。若假設(shè)不成立：可以嘗試變換協(xié)變量或因變量或使用能對每個斜率獨(dú)立解釋的模型或使用不需要假設(shè)回歸斜率同質(zhì)性的非參數(shù)ancova方法。結(jié)果可視化檢驗(yàn)斜率同質(zhì)性的可視化方法install.packages("HH")library(HH)ancova(weight~gesttime+dose,data=litter)雙因素方差分析例子：60只豚鼠，分別采用兩種喂食方法（橙汁及維C），各種喂食方法中抗壞血酸含量有三種水平（0.5mg\1mg\2mg）,每種處理方式分配給10只豚鼠。因變量為牙齒長度。主效應(yīng)和交互效應(yīng)都非常顯著意味著什么？即喂食方式、喂食劑量、以及二者交互都影響方差一致性。是嗎？可視化處理：帶置信區(qū)間的均值以上只顯示了交互效應(yīng)，未顯示主效應(yīng)。推薦下圖方法：library(HH)interaction2wt(len~supp*dose)模型的假設(shè)檢驗(yàn)及均值比較及上節(jié)類似。重復(fù)測量方差分析受試者不止被測量一次。一般情況含一個組間因子，一個組內(nèi)因子。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="Carbondioxideuptakerate(umol/m^2sec)")二氧化碳濃度及植物類型對二氧化碳吸收的交互效應(yīng)解釋：ChilledQuebec比MississippiPlants二氧化碳吸收能力強(qiáng)，隨著濃度增加，這種差異越來越明顯。多元方差分析當(dāng)因變量不止一個時，可用多元方差分析。單因素多元方差分析例子：研究鼓舞中卡路里、脂肪、糖是否會因?yàn)閮Υ嫖恢茫?代表底層貨架，2中層，3高層）不同而發(fā)生變化?？防铩⒅?、糖為因變量，儲存位置為自變量。上面僅說了不同貨架上含卡路里脂肪糖均是不同的。可以補(bǔ)充個均值比較（如TukeyHSD）來判斷對于每個因變量，哪個貨架及其他貨架是不同的。評估假設(shè)檢驗(yàn)單因素多元方差分析假設(shè)：多元正態(tài)性向量x及均值u.x及均值的馬氏距離平方服從自由度為p的卡方分布該例子中。P=3QQ圖展示卡方分布的分位數(shù)，橫縱坐標(biāo)分別是樣本量及馬氏距離平方值。方差—協(xié)方差矩陣同質(zhì)性：各組的協(xié)方差矩陣相同。通?？捎肂ox’sM來評估該假設(shè)。但R中么有。對多元正態(tài)性敏感，很容易拒絕該假設(shè)。檢驗(yàn)多元離群點(diǎn)library(mvoutlier)outliers<-aq.plot(y)outliers穩(wěn)健多元方差分析如果多元正態(tài)性或者方差—協(xié)方差均值假設(shè)都不滿足，又擔(dān)心多元離群點(diǎn)，可以用穩(wěn)健或者非參數(shù)版的MANOVA檢驗(yàn)。不知道具體用啥模型用回歸做ANOVA沒看十、功效分析主要內(nèi)容：判斷所需樣本量計(jì)算效應(yīng)值評價統(tǒng)計(jì)功效樣本大?。簩τ谘芯坑蠿個受試者，值得做嗎？到底需要多少個受試者顯著性：H0為真，拒絕H0的概率（1型錯誤）功效：H0假設(shè)為假時，將H0拒絕的概率：（1-11型錯誤）效應(yīng)值：依據(jù)所用假設(shè)檢驗(yàn)不同，如：t檢驗(yàn)，效應(yīng)值為d=假設(shè)檢驗(yàn)概述功效分析T檢驗(yàn)方差分析相關(guān)性線性模型比例檢驗(yàn)卡方檢驗(yàn)效應(yīng)值功效分析圖形十一、中級繪圖1散點(diǎn)圖散點(diǎn)圖矩陣高密度散點(diǎn)圖（重合點(diǎn)很多）用封箱、顏色、透明度來指明任意點(diǎn)上重疊點(diǎn)的數(shù)目封箱圖密度通過顏色標(biāo)識三維散點(diǎn)圖氣泡圖一般通過三維圖展示三個定量變量之間關(guān)系，還有一種思路：先創(chuàng)建一個二維散點(diǎn)圖，然后用點(diǎn)的大小來代表第三個變量的值，這便是氣泡圖。Plot()調(diào)用即添加新圖Line()在原圖形上添加信息繪圖參數(shù)亂了，重啟就恢復(fù)默認(rèn)了2折線圖#代碼清單11-4展示五種橘樹隨時間推移的生長狀況的折線圖Orange$Tree<-as.numeric(Orange$Tree)ntrees<-max(Orange$Tree)xrange<-range(Orange$age)yrange<-range(Orange$circumference)plot(xrange,yrange,type="n",xlab="Age(days)",ylab="Circumference(mm)")colors<-rainbow(ntrees)linetype<-c(1:ntrees)plotchar<-seq(18,18+ntrees,1)for(iin1:ntrees){tree<-subset(Orange,Tree==i)lines(tree$age,tree$circumference,type="b",lwd=2,lty=linetype[i],col=colors[i],pch=plotchar[i])title("TreeGrowth","exampleoflineplot")legend(xrange[1],yrange[2],1:ntrees,cex=0.8,col=colors,pch=plotchar,lty=linetype,title="Tree"比較五種橘樹生長狀況的折現(xiàn)圖3相關(guān)圖順