第四章R編程實(shí)踐

1、第四章 R 編程實(shí)踐自此，我們已經(jīng)對(duì)R 的功能有個(gè)全面的了解。一些數(shù)據(jù)分析都是在R 對(duì)話窗口( R Console )中進(jìn)行的，但對(duì)于復(fù)雜的數(shù)據(jù)分析顯然是不行方便的。下面將從統(tǒng)計(jì)語(yǔ)言和編程角度來(lái)說(shuō)明，我們將會(huì)了解一些在 R編程實(shí)踐中采用的簡(jiǎn)單的概念。一、 循環(huán)和向量化相比下拉菜單式的程序(R下開發(fā)的統(tǒng)計(jì)包R Commander, R的一個(gè)優(yōu)勢(shì)在于它可以把一系列連續(xù)的操作簡(jiǎn)單的程序化。這一點(diǎn)上和所有其他計(jì)算機(jī)編程語(yǔ)言是一致的，但R 有一些特性使得非專業(yè)人士也可以很簡(jiǎn)單的編寫程序。和其他編程語(yǔ)言一樣，R有一些和C語(yǔ)言類似的控制結(jié)構(gòu)。(一)控制結(jié)構(gòu)1. 條件語(yǔ)句：條件語(yǔ)句常用于避免除零或負(fù)數(shù)等數(shù)學(xué)問(wèn)

2、題。它有兩種形式：1)if ( 條件 ) 表達(dá)式 1 else 表達(dá)式 22)ifelse( 條件 ,yes,no)例如：if (x=0) sqrt(x) else NAifelse (x=0,sqrt(x),NA)2. 循環(huán)：有三種形式：1)for ( 變量 in 向量 ) 表達(dá)式2)while( 條件 ) 表達(dá)式3)repeat 表達(dá)式例如：假定我們有一個(gè)向量x ，對(duì)于向量x 中值為 b 的元素，我們把0 賦給另外一個(gè)等長(zhǎng)向量y 的對(duì)應(yīng)元素，否則賦1。我們首先創(chuàng)建一個(gè)與向量x 等長(zhǎng)的向量y ：y - numeric(length(x)for (i in 1:length(x) if (xi

3、 = b) yi - 0 else yi - 1幾個(gè)指令可以放在一個(gè)大括弧里面：for (i in 1:length(x) yi - 0x - rnorm(10,-5,0.1)y - rnorm(10,5,2)X - cbind(x,y) # 矩陣X的列名保持x和y for (i in 1:5/2) zi minimum) count-0while (count10) print(count)count=count+1#如果設(shè)置不正確while 循環(huán)可以能導(dǎo)致無(wú)限循環(huán)，要小心#computeEstimate()x0-1tol-1e-5repeatx1-rnorm(1,1,2)if (abs(x

4、1-x0)tol)breakelsex0-x1#next,returnfor(i in 1:100)if(i=20)#skipnext#do somthing here3. 分支控制語(yǔ)句switch(statement,list)switch(2, banana , apple , other )num-10mode z z for (i in 1:length(z) zi yx= b yx!= b x y X apply(X,2,mean)lapply() 可以用于一個(gè)列表對(duì)象：它的語(yǔ)法類似apply 并且返回一個(gè)列表對(duì)象。 forms lapply(forms,lm)sapply() 是函

5、數(shù) lapply() 一個(gè)更為靈活的變種，它可以接受向量或者矩陣作為主要參數(shù)，并且返回形式更為友好的結(jié)果，常常是表格方式。還有sweep() 、 tapply() 函數(shù)等。例如：attach(mtcars)tapply(mpg,cy1,mean)col.mean=apply(mtcars,2,mean)sweep(mtcars,2,col.mean,-)表示從數(shù)據(jù)框mtcars 中按列計(jì)算col.mean ，并從 mtcars 中減去。2、 用 R 寫程序一般情況下，一個(gè) R程序以ASCII格式保存，擴(kuò)展名為.R。如果一個(gè)工作要重復(fù)好多次，用 R 程序是一個(gè)不錯(cuò)的選擇。在我們的第一個(gè)例子中，我

6、們想對(duì)三個(gè)不同種屬的鳥繪制一樣的圖，而 且數(shù)據(jù)在三個(gè)不同的文件中。我們將一步一步的演示看R用不同的方式去完成這個(gè)簡(jiǎn)單的問(wèn)題。首先，我們憑直覺連續(xù)鍵入一系列命令，而且預(yù)先分割圖形設(shè)備。layout(matrix(1:3,3,1) # 分割圖形界面data - read.table(Swal.dat) # 讀入數(shù)據(jù)plot(data$V1,data$V2,type=l)title(swallow) # 增加標(biāo)題data - read.table(Wren.dat)plot(data$V1,data$V2,type=l)title(wren)data - read.table(Dunn.dat)pl

7、ot(data$V1,data$V2,type=l)title(dunnock)字符#用于在程序中添加注釋行：R會(huì)自動(dòng)跳過(guò)注釋行。第一個(gè)程序的問(wèn)題是在我們加入另外一個(gè)物種數(shù)據(jù)時(shí)，它過(guò)長(zhǎng)。此外，一些命令多次執(zhí)行， 因 此它們可以放在一起，在執(zhí)行的時(shí)候僅僅修改一些參數(shù)。這里的策略是把參數(shù)放到一個(gè)字符型的向 量中去，然后用下標(biāo)去訪問(wèn)這些不同的值。layout(matrix(1:3,3,1) # 分割圖形界面species - c(swallow,wren,dunnock)file - c(Swal.dat,Wren.dat,Dunn.dat)for(i in 1:length(species) da

8、ta - read.table(filei) #讀入數(shù)據(jù)plot(data$V1,data$V2,type=l) title(speciesi) # 增加標(biāo)題 注意代碼 read.table() 里面的參數(shù)filei 上面沒有雙引號(hào)， 因?yàn)檫@個(gè)參數(shù)是字符型?，F(xiàn)在的代碼比較緊湊。 它比較容易加入新的物種， 因?yàn)樵O(shè)置物種名字和數(shù)據(jù)文件的向量都程序 的前端。如果擴(kuò)展名為.dat的數(shù)據(jù)文件在 R的工作目錄下面，程序可以正常運(yùn)行，否則用戶要設(shè)置工作目錄，或者設(shè)置絕對(duì)路徑( 例如： file source(Mybirds.R)和所有以文件作為輸入對(duì)象的函數(shù)一樣，如果該文件不在當(dāng)前工作目錄下面，用戶需要提供

9、該文件的絕對(duì)路徑。3、 編寫你自己的函數(shù)(一)編寫自己的函數(shù)1. 函數(shù)的基本操作大多數(shù) R 的工作是通過(guò)函數(shù)來(lái)實(shí)現(xiàn)，而且這些函數(shù)的輸入?yún)?shù)都放在一個(gè)括弧里面。用戶可以編寫他們自己的函數(shù)，并且這些函數(shù)和R里面的其他函數(shù)有一樣的特性。stdev=function(x)# 計(jì)算標(biāo)準(zhǔn)差 sqrt(var(x)stdevz=1:3stdev(z)函數(shù)的返回值可以是任何R對(duì)象，甚至可以是另一個(gè)函數(shù)?？梢酝ㄟ^(guò)顯式地調(diào)用return(),把一個(gè)值返回給主調(diào)函數(shù)。如果不用這條語(yǔ)句，默認(rèn)會(huì)將最后執(zhí)行的語(yǔ)句的值作為返回值。 例如：oddcount-function(x)k-0#k 的初始值 for(n in x)i

10、f(n%2=1)k=k+1#%是模 return(k)這個(gè)函數(shù)返回參數(shù)x 中奇數(shù)的個(gè)數(shù)。如果函數(shù)的輸出包括多個(gè)變量，可以把他們集中在一個(gè)列表中，以列表為參數(shù)調(diào)用函數(shù)，讓這個(gè)函數(shù)返回列表。創(chuàng)建函數(shù)之后可以通過(guò)函數(shù)formals() 和 body() 來(lái)查看其參數(shù)。formals(oddcount)body(oddcount)可以使用 page() 一個(gè)函數(shù)的代碼，確定它的功能和用法。page(abline)函數(shù)內(nèi)變量的賦值只在函數(shù)內(nèi)有效。編寫自己的函數(shù)可以讓你有效的，靈活的與合理的使用R。我們?cè)俅问褂们懊孀x數(shù)據(jù)并且畫圖的例子。如果我們想在其他情況下進(jìn)行這樣的操作，寫一個(gè)函數(shù)是一個(gè)不錯(cuò)的想法：my

11、fun - function(S,F) data - read.table(F)plot(data$V1,data$V2,type=l)title(S) 執(zhí)行時(shí)，這個(gè)函數(shù)必須載入內(nèi)存。當(dāng)然，這有多種方式實(shí)現(xiàn)。和所有其他命令一樣，函數(shù)的各行可以直接通過(guò)鍵盤鍵入， 或者從一個(gè)文本編輯器里面拷貝粘貼。 如果函數(shù)保存在一個(gè)文本文件中，可以命令 source() 載入。如果用戶期望一些函數(shù)在 R 啟動(dòng)時(shí)就被載入，可以把它們保存在工作目錄下面的文件.RData中。另外一種方式是，配置文件.Rpro le或Rpro - le( 詳見?Startup fordetails) 。最后，還可以創(chuàng)建一個(gè)包，但是這里

12、不想多討論。一旦函數(shù)載入后，我們就可以鍵入一條命令以讀入數(shù)據(jù)和畫出我們想要的圖，如myfun(swallow,Swal.dat) 。因此，現(xiàn)在我們的程序有第三個(gè)實(shí)現(xiàn)的版本了：layout(matrix(1:3,3,1)myfun(swallow,Swal.dat)myfun(wren,Wrenn.dat)myfun(dunnock,Dunn.dat)我們還可以用 sapply() 實(shí)現(xiàn)程序的第四個(gè)版本：layout(matrix(1:3,3,1)species - c(swallow,wren,dunnock)file foo x foo()1 1x 不是為了在函數(shù)foo() 里面創(chuàng)建對(duì)象，

13、因此 R 將會(huì)在封裝環(huán)境中搜索是否有個(gè)名為 x 的對(duì)象，和打印它的值( 否則一條錯(cuò)誤信息將會(huì)顯示，執(zhí)行中斷) 。如果 x 是我們定義的函數(shù)中一個(gè)對(duì)象的名字，全局環(huán)境中變量x 值將會(huì)被采用。 x foo2 - function() x foo2()1 2 x1 1此時(shí)， print() 使用在它所在的環(huán)境中定義的x ，即環(huán)境foo2 中的 x 。前面提及的 封裝 是關(guān)鍵所在。在前面兩個(gè)演示函數(shù)中，有兩個(gè)環(huán)境：全局環(huán)境和函數(shù)foo 或foo2 的局部環(huán)境。 如果有三個(gè)或者更多的嵌套環(huán)境，對(duì)象搜索將逐層搜索直到全局環(huán)境。rm(list=ls()ls()ls.str() 可以獲得更多信息。有兩種方式指

14、定一個(gè)函數(shù)的參數(shù)： 通過(guò)它們的定義時(shí)的位置或者名字( 又稱為標(biāo)簽參數(shù)) 。 例如，假定一個(gè)函數(shù)有三個(gè)參數(shù)：foo - function(arg1,arg2,arg3) foo()在執(zhí)行時(shí)可以不用名字argl,，如果相應(yīng)的參數(shù)對(duì)象放在相應(yīng)的位置上，如：foo(x,y,z) 。但是，如果使用了參數(shù)的名字，位置信息將會(huì)失效，如 foo(arg3 = z,arg2 = y,arg1 = x) 。 R 函數(shù)的另外一個(gè)特性是函數(shù)可能采用定義時(shí)的默認(rèn)設(shè)置。例如：foo - function(arg1,arg2 = 5,arg3 = FALSE) 命令 foo(x) ， foo(x,5,FALSE) 和 fo

15、o(x,arg3 = FALSE) 將會(huì)得到一樣的結(jié)果。使用一個(gè)函數(shù)的默認(rèn)設(shè)置非常有用，特別在使用標(biāo)簽參數(shù)的時(shí)候，如 foo(x,arg3 = TRUE) 僅僅改變一個(gè)默認(rèn)設(shè) 置。環(huán)境層次中，某一層次的代碼對(duì)它上級(jí)層次中所有變量至少有讀的權(quán)限。但是，通過(guò)標(biāo)準(zhǔn)的運(yùn) 算符 -直接對(duì)上級(jí)層次的變量進(jìn)行寫操作是不可行的。#觀察 w 的值是否有變化w-12f-function(y)d-8w-w+1y-y-2cat(w=,w,n)h-function() return(d*(w+y) return(h() f(2) cat(w=,w,n) 如果你希望對(duì)一個(gè)全局變量，或者更一般情況下，對(duì)當(dāng)前層次的上級(jí)層次中

16、任意變量進(jìn)行寫操 作，可以在當(dāng)前層次使用超賦值運(yùn)算符- ，或者函數(shù)assign() 。(二)函數(shù)的調(diào)用函數(shù)也是對(duì)象，也可以給它們賦值，把它們用作其他函數(shù)的參數(shù)，如下所示：#函數(shù)的調(diào)用f1-function(a,b) return(a+b)f2-function(a,b) return(a-b)f-f1f(3,2)f-f2f(3,2)g-function(h,a,b) h(a,b)g(f1,3,2)g(f2,3,2)#在同一幅畫上繪制若干個(gè)函數(shù)的圖形g1-function(x) return(5*sin(x)+5)g2-function(x) return(sqrt(xA2+1)g3-funct

17、ion(x) return(2*x-1)plot(c(0,6.5),c(-5,15)# 準(zhǔn)備繪圖，確定 x和y的區(qū)域for(f in c(g1,g2,g3)plot(f,0,6.5,add=T)#在原有的圖形上繪圖(三)函數(shù)編輯R提供了一個(gè)edit函數(shù)，讓使用者可以在一個(gè)特定的窗口中更改自己的函數(shù)；fix使用預(yù)設(shè)的編輯器。fix(stdev)edit(file=for.r)我們來(lái)看另外一個(gè)例子。盡管這個(gè)例子不是純的統(tǒng)計(jì)學(xué)例子，但是它很好地展示了R語(yǔ)言的靈活性。假定我們想研究一個(gè)非線性模型的行為：這個(gè)模型( Ricker模型)的定義如下：Nt 1 Nt exp r 1 N K這個(gè)模型廣泛地用于種

18、群動(dòng)態(tài)變化( population dynamics )的研究里面，特別是魚類的種群變 化。我們想用一個(gè)函數(shù)去模擬這個(gè)模型關(guān)于增長(zhǎng)率r和初始群體大小 N的變化情況(承載能力K常常設(shè)定為1且以這個(gè)值作為默認(rèn)值)的影響；結(jié)果將以種群大小相對(duì)時(shí)間的圖表示。我們還將設(shè)定一個(gè)可選項(xiàng)允許用戶只顯示最后若干步中種群大小(默認(rèn)所有結(jié)果都會(huì)被繪制出來(lái))。下面的函數(shù)就是做Ricker模型的數(shù)值模擬的。ricker - function(nzero,r,K=1,time=100,from=0,to=time) N - numeric(time+1)N1 - nzerofor (i in 1:time) Ni+1 -

19、 Ni*exp(r*(1 - Ni/K)Time layout(matrix(1:3,3,1) ricker(0.01,0.1);title(r = 0.1) ricker(0.05,0.2);title(r = 0.2) ricker(0.1,0.3);title(r = 0.3)四、防錯(cuò)、除錯(cuò)與調(diào)試1 .防錯(cuò)與除錯(cuò)R具備在函數(shù)中加入防錯(cuò)功能的應(yīng)用函數(shù)：stop()和warning。stop()會(huì)傳回錯(cuò)誤信息并終止執(zhí)行函數(shù)；warning。則在傳回錯(cuò)誤訊息的同時(shí)繼續(xù)執(zhí)行函數(shù)。stdev=function(x) #計(jì)算標(biāo)準(zhǔn)差if (!is.vector(x)stop(向量才能計(jì)算標(biāo)準(zhǔn)差)7el

20、se sqrt(var(x)z=3z=as.matrix(z)stdev(z)# 利用輸出cat 語(yǔ)句進(jìn)行確認(rèn)g=function(x,y)xA2-y+1z-0x-yA2+3*g(z,2)cat(x=,x,n)#利用print 語(yǔ)句進(jìn)行調(diào)試x-yA2+3*g(z,2)cat(x=,x,n)w-28q0)u-1print(the if was done)elsev-10print(the else was done)2. 調(diào)試R 的核心調(diào)試工具由瀏覽器(browser )構(gòu)成，它可以讓你逐行運(yùn)行代碼， 并在運(yùn)行過(guò)程中檢查?？梢杂糜袃?nèi)置除錯(cuò)函數(shù)debug() 或 browser()函數(shù)來(lái)打開這一瀏

21、覽器。R 的調(diào)試工具是針對(duì)單個(gè)函數(shù)的 ，如果你認(rèn)為函數(shù) f() 中有漏洞，就可以調(diào)用 debug(f) 來(lái)設(shè)置 函數(shù)的調(diào)試狀態(tài)。調(diào)用undebug(f) 將取消函數(shù)的調(diào)試狀態(tài)。debug(stdev)stdev(z)如果你只想對(duì)f() 進(jìn)行一次調(diào)試會(huì)話，可以使用 debugonce(f) 函數(shù)，會(huì)在f() 第一次執(zhí)行時(shí)對(duì)其設(shè)置調(diào)試狀態(tài)。你果你認(rèn)為漏洞離函數(shù)的開始部分較遠(yuǎn)，可以在函數(shù)的某一行加入 browser() 語(yǔ)句，那么瀏覽器只在程序執(zhí)行到這一行時(shí)被打開。#設(shè)置斷點(diǎn)f=function(x)xA3-2*x-1root=function(x,y)if(f(x)*f(y)=0)browser(

22、)x1-x;x20)stop( 區(qū)間端點(diǎn)的 f() 值之正負(fù)號(hào)不相反)elseif(f(x)0) x1-y;x2-xrepeatx00)x2-x0elsex1-x0if (abs(x1-x2)為變更為Browserd ( d 表示函數(shù)調(diào)用鏈的深度) 。在命令提示符之后你可以輸入以下命令：n (代表 next ) ：執(zhí)行下一行然后暫停，直接鍵入 Enter 也一樣。c (代表continue ):與n類似，但是如果當(dāng)前在某個(gè)循環(huán)中，將會(huì)執(zhí)行剩余的循環(huán)退出循環(huán)才會(huì)暫停。如果你不在某個(gè)循環(huán)中，那么下一次暫停前函數(shù)的剩余部分將被執(zhí)行。任意的R命令：在瀏覽器中，你依然處在R的交互模式，因此可以任意查看變量的取值。where ：會(huì)輸出一份棧跟蹤路徑，即顯示到達(dá)當(dāng)前