R語言生物統(tǒng)計

1、R語言概述 R是屬于GNU系統(tǒng)的一個自由、免費、源代碼開放的軟件，它是一個用于統(tǒng)計計算和統(tǒng)計制圖的優(yōu)秀工具。R是S語言的一種實現(xiàn)。S語言是由AT&T貝爾實驗室開發(fā)的一種用來進行數(shù)據(jù)探索、統(tǒng)計分析、作圖的解釋型語言。最初S語言的實現(xiàn)版本主要是S-PLUS。S-PLUS是一個商業(yè)軟件，它基于S語言，并由MathSoft公司的統(tǒng)計科學部進一步完善。后來Auckland大學的Robert Gentleman 和 Ross Ihaka 及其他志愿人員開發(fā)了一個R系統(tǒng)。R的使用與S-PLUS有很多類似之處，兩個軟件有一定的兼容性。S-PLUS的使用手冊，只要經(jīng)過不多的修改就能成為R的使用手冊。所以

2、有人說：R，是S-PLUS的一個“克隆”。但是請不要忘了：R is free。 R是一套完整的數(shù)據(jù)處理、計算和制圖軟件系統(tǒng)。其功能包括：數(shù)據(jù)存儲和處理系統(tǒng)；數(shù)組運算工具（其向量、矩陣運算方面功能尤其強大）；完整連貫的統(tǒng)計分析工具；優(yōu)秀的統(tǒng)計制圖功能；簡便而強大的編程語言：可操縱數(shù)據(jù)的輸入和輸入，可實現(xiàn)分支、循環(huán)，用戶可自定義功能。 與其說R是一種統(tǒng)計軟件，還不如說R是一種數(shù)學計算的環(huán)境，因為R并不是僅僅提供若干統(tǒng)計程序、使用者只需指定數(shù)據(jù)庫和若干參數(shù)便可進行一個統(tǒng)計分析。R的思想是：它可以提供一些集成的統(tǒng)計工具，但更大量的是它提供各種數(shù)學計算、統(tǒng)計計算的函數(shù)，從而使使用者能靈活機動的進行數(shù)據(jù)分

3、析，甚至創(chuàng)造出符合需要的新的統(tǒng)計計算方法。 R是一個免費的自由軟件，它有UNIX、LINUX、MacOS和WINDOWS版本，都是可以免費下載和使用的。R的主要網(wǎng)站是。在那兒可以下載到R的安裝程序、各種外掛程序和文檔。在R的安裝程序中只包含了8個基礎模塊，其他外在模塊可以通過CRAN獲得（）。 以下簡述R FOR WINDOWS的安裝和使用： 在網(wǎng)址： 下可以找到R的各個版本的安裝程序和源代碼。點擊進入：Windows (95 and later)，再點

4、擊：base，下載SetupR.exe，約18兆，此便是R FOR WINDOWS的安裝程序。雙擊SetupR.exe，按照提示一步步安裝即可。 安裝完成后，程序會創(chuàng)建R程序組并在桌面上創(chuàng)建R主程序的快捷方式（也可以在安裝過程中選擇不要創(chuàng)建）。通過快捷方式運行R，便可調(diào)出R的主窗口（如下圖 1-1）。 類似于許多以編程方式為主要工作方式的軟件，R的界面簡單而樸素，只有不多的幾個菜單和快捷按鈕。快捷按鈕下面的窗口便是命令輸入窗口，它也是部分運算結(jié)果的輸出窗口，有些運算結(jié)果則會輸出在新建的窗口中。 主窗口上方的一些文字是剛運行R時出現(xiàn)的一些說明和指引。文字下的：> 符號便是R的命令提示符，在

5、其后可輸出命令；>后的矩形是光標。R一般是采用交互方式工作的，在命令提示符后輸入命令，回車后便會輸出結(jié)果。第一節(jié) 數(shù)據(jù)的輸入注：在R中，大小寫是敏感的，即大小寫是不同的。  一、變量數(shù)據(jù)直接輸入  1、數(shù)值變量  如果有一個變量，變量名為x，其數(shù)據(jù)如下：12，15，46，23，15。命令語句如下：  x<-c(12,15,46,23,15)  或   c(12,15,46,23,15)->x  其中：x為變量名；<-與->為賦值符；c( )為向量

6、建立函數(shù)，表示把括號中的數(shù)據(jù)建立為一個向量。  以上命令語句建立了一個數(shù)值變量x。 2、字符變量  字符變量的建立與數(shù)值變量一致。字符串使用引號（單、雙均可），如： y<-c(“er”,”sdf”,”eir”,”jk”,”dim”)   或   c(“er”,”sdf”,”eir”, ”jk”,”dim”)->y 將建立字符變量y。 3、邏輯變量  邏輯變量中的元素是：TRUE(或簡寫為T)、FALSE(F)、NA(表示缺省)。請注意必須都是大

7、寫。  邏輯變量可以直接輸入，如： z<-c(TRUE,FALSE,F,T,NA) 可建立邏輯變量z。  1可以表示T，0可以表示F。所以可以直接寫成： z<-c(1，0，0，1,NA)  在計算時，T和F也是當作1和0來用的。 4、其他輸入方式  1）利用“：”      “：”可以產(chǎn)生連續(xù)的整數(shù)序列，如     > x<-1:10   > x&

8、#160;   1  1  2  3  4  5  6  7  8  9 102）利用seq函數(shù)  產(chǎn)生有規(guī)律的數(shù)列的更一般方法是采用seq函數(shù)，它的一般用法是：seq（下界，上界，間距），如： > x=seq(1,10,0.5)  > x    1  1.0

9、0; 1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0  8.5  9.0  9.5 10.0  3）利用rep函數(shù)  rep函數(shù)重復第一個變量若干次&

10、#160;> y<-rep(1,5)  > y 1 1 1 1 1 1  二、建立數(shù)據(jù)集  在上面，我們建立起了三個變量：x，y，z。如果要建立起一個三變量的數(shù)據(jù)集，即含有x、y、z的二維表，則可使用data.frame()函數(shù)。當然，在這種情況下，要求變量的長度是一樣的。  dd<-data.frame(x,y,z) 便可建立數(shù)據(jù)集dd。 edit(dd)  可調(diào)用R中的數(shù)據(jù)編輯器顯示、編輯數(shù)據(jù)集dd。如下圖。

11、     如果eidt( )編輯的是一維的向量，如：edit(x)，那么它便會調(diào)用記事本來編輯，而不會調(diào)用數(shù)據(jù)編輯器。如下圖。     三、從其他文件中讀取數(shù)據(jù)  在應用統(tǒng)計學中，數(shù)據(jù)量一般是比較大的，變量也很多。如果用上述方法來建立數(shù)據(jù)集，好像辛苦了一些。上述方法適用于少量數(shù)據(jù)、少量變量的分析。對于大量數(shù)據(jù)和變量，一般應在其他軟件中輸好，再讀R中處理。 1、 讀入輸好的純文本文件  若數(shù)據(jù)已經(jīng)輸入一個純文本文件“c:/test.txt”。其格式如下：     &#

12、160;Price    Floor     Area   Rooms     Age  Cent.heat 01   52.00    111.0      830     5   

13、60;   6.2      no 02   54.75    128.0      710     5       7.5      no 03   

14、;57.50    101.0     1000     5       4.2    no 04   57.50    131.0      690     

15、6       8.8      no 05   59.75     93.0      900     5       1.9      yes&#

16、160; 其中第一行為變量名，第一列為記錄序號。 則利用read.table( )函數(shù)可讀入數(shù)據(jù)，如： read.table(“c:/test.txt”)->test  便把數(shù)據(jù)集讀入，并命名為test。 注：文件的后綴不必一定要.txt，關鍵文件要為純文本，里面不能有特殊格式符。 如果數(shù)據(jù)文件中沒有第一列記錄序號，如：Price Floor Area Rooms Age Cent.heat 52.00 111.0 830 5 6.2 no 54.75 128.0 710 5 7.5 no 57.50 101.0 1000 5

17、4.2 no 57.50 131.0 690 6 8.8 no 59.75 93.0 900 5 1.9 yes 則命令語句為： read.table(“c:/test.txt”, header=TRUE)->test 系統(tǒng)根據(jù)每個變量第一個值的類型，自動識別變量類型，如以上數(shù)據(jù)集中，除最后一個“Cent.hear”是字符變量，其他均為數(shù)值變量。 2、 讀入其他格式的數(shù)據(jù)庫 要讀入其他格式數(shù)據(jù)庫，必須先安裝“foreign”模塊。它不屬于R的8個內(nèi)在模塊，需要在使用前安裝。 安裝方法很簡便，只需鍵入命令： library(foreign) 即可。 四、向量運算 1、加減乘除 對向量進行加

18、減乘除運算，其含義是對向量的每一個元素進行運算： 例： > x<-c(1,4,7,8,10) > y<-x+2 > y 1 3 6 9 10 12 向量之間的運算，則是對應元素之間的運算。 例： > x<-c(1,4,7,8,10) > y<-x+2 > z<-x+y > z 1 4 10 16 18 222、 添加數(shù)據(jù)與替換數(shù)據(jù) 添加函數(shù)的句法是：append(向量，添加部分，after=起點)。 例： > x<-seq(1,15,2) > x 1 1 3 5 7 9 11 13 15 > app

19、end(x,20:30,after=4) 1 1 3 5 7 20 21 22 23 24 25 26 27 28 29 30 9 11 13 15 如果after參數(shù)缺省，則默認在向量的最后插入。 > x<-seq(1,15,2) > append(x,20:30) 1 1 3 5 7 9 11 13 15 20 21 22 23 24 25 26 27 28 29 30 替換數(shù)據(jù)的句法是：replace（向量，替換位置，替換值）。 例： > x<-seq(1,15,2) > replace(x,c(1,4,7),-1) 1 -1 3 5 -1 9 11

20、-1 15 第二節(jié) 統(tǒng)計描述一、 數(shù)值變量的統(tǒng)計描述 在R中沒有一個直接的統(tǒng)計描述命令，可以把常見的指標都一起算出來。如概述中所說，R往往是提供一些比較基礎的統(tǒng)計命令，需要使用者自己編程組合。 1、和：sum（x），積：prod（） 例： > x<-c(12,14,16,18) （建立一個數(shù)據(jù)向量x） > sum(x) （求x的和） 1 60。 > prod(x) 1 48384 2、算術均數(shù)：mean(x) 例： > x<-c(12,14,16,18) （建立一個數(shù)據(jù)向量x） > mean(x) （求x的算術平均數(shù)） 將輸出x的均數(shù)15。 但是，算術

21、平均數(shù)有一個缺點就是對outlier非常敏感，如果少數(shù)幾個值非常大，那么整個算術平均數(shù)也隨之增加。因此，為了避免這種情況，可以使用trimmed mean也就是說，先對最大和最小的數(shù)剔除了（一般會按照百分比剔除，比如剔除最大和最小的a的個案剔除出去），然后再計算算術平均值。mean命令也可以計算trimmd mean，只要增加一個trim參數(shù)就可以了。 > c(1,2,3,4,5,6,7,8,9,12,1,31,4,15,16,17)->z1 > mean(z1) 1 8.8125 > mean(z1,trim=0.1) （剔除最大和最小10%的數(shù)據(jù)） 1 7.7857

22、14 有時候，計算mean時的權(quán)重希望不是一樣的，那么可以使用weighted.mean命令計算。比如在評委評分的時候，可能資深評委的評委更重要，給予他們更大的權(quán)重就是這種情況。比如有五位評委，他們的某一次評分分別是2，3，4，5，4，而他們的權(quán)重是0.1，0.2，0.2，0.3，0.6那么可以使用以下的命令求出加權(quán)平均數(shù)。 > weighted.mean(x=c(2,3,4,5,4),w=c(0.1,0.2,0.2,0.3,0.6) 13.9285713、中位數(shù)：median(x) 例： > y<-c(11,13,14,12,17,10) > median(y) 將輸

23、出中位數(shù)12.5 4、極差：max(x)-min(x) 例： > y<-c(11,13,14,12,17,10) > y=max(y)-min(y) 將輸出結(jié)果7。 還有一個命令是range命令，他同時返回最大值和最小值。 > y<-c(11,13,14,12,17,10) > range(y) 1 10 17 5、樣本方差：var(x) 例： > x<-c(9,8,7,10,12,10,11,14,8,9) > var(x) 將輸出結(jié)果4.4。 6、樣本標準差：sd(x) 例： > x<-c(9,8,7,10,12,10,11

24、,14,8,9) > sd(x) 將輸出結(jié)果2.097618。 請注意都是“樣本”。 7、變異系數(shù)：(sd(x)/mean(x)*100 8、四分位差：IQR（x） > x<-c(9,8,7,10,12,10,11,14,8,9) > IQR(x) 輸出結(jié)果 2.5。 注：四分位差表示的，是75分位值和25分位值的差 通過以上命令和程序，可以求出想要的描述指標。 二、正態(tài)性檢驗與方差齊次性檢驗 1正態(tài)性檢驗 可采用Shapiro-Wilk檢驗。程序如下： > x<-c(9,8,7,10,12,10,11,14,8,9) > shapiro.test(x

25、) 輸出結(jié)果如下： Shapiro-Wilk normality testdata: x W = 0.9513, p-value = 0.68432方差齊次性分析 例： > y<-c(10,12,14,15,17,29) > x<-c(9,8,7,10,12,10,11,14,8,9) > var.test(x,y) 結(jié)果如下： F test to compare two variances data: x and y F = 0.097, num df = 9, denom df = 5, p-value = 0.003206 alternative hypot

26、hesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.01451680 0.43493188 sample estimates: ratio of variances 0.0969875 第三節(jié) t檢驗 一、 樣本均數(shù)與總體均數(shù)的假設檢驗 例：隨機抽測了10只兔的直腸溫度，其數(shù)據(jù)為：38.7、39.0、38.9、39.6、39.1、39.8、38.5、39.7、39.2、38.4（），已知該品種兔直腸溫度的總體平均數(shù)0µ=39.5（），試檢驗該樣本平均溫度與0µ是

27、否存在顯著差異？ R程序如下： > x<-c(38.7,39,38.9,39.6,39.1,39.8,38.5,39.7,39.2,38.4) （建立變量x） > t.test(x,mu=39.5,alternative="two.sided",conf.level=0.95) 注：以上程序中，x為數(shù)據(jù)變量；mu=39.5為總體均數(shù)；alternative=“two.sided”表示行雙側(cè)檢驗，如果alternative=“l(fā)ess”，為1<2的單側(cè)檢驗，反之，alternative=“greater”，為1>2的單側(cè)檢驗；conf.level

28、=0.95指檢驗水準為0.05。 輸出結(jié)果如下： One Sample t-test data: x t = -2.641, df = 9, p-value = 0.02687 alternative hypothesis: true mean is not equal to 39.5 95 percent confidence interval: 38.73882 39.44118 sample estimates: mean of x 39.09 輸出假設檢驗結(jié)果、備擇假設、總體均數(shù)95可信區(qū)間和樣本均數(shù)。從結(jié)果中可以看出，t=-2.641，自由度為9，樣本均數(shù)95%的可信區(qū)間為38.73

29、882，39.44118，假設檢驗的概率為0.02687，因此，0.01<P<0.05，備擇假設成立，兩者不相等，差異顯著，樣本均數(shù)為39.09。 二、配對資料T檢驗 例：11只60日齡的雄鼠在x射線照射前后之體重數(shù)據(jù)見下表（單位：g）：檢驗雄鼠在照射x射線前后體重差異是否顯著？ 編 號 1 2 3 4 5 6 7 8 9 10 11 照射前 25.7 24.4 21.1 25.2 26.4 23.8 21.5 22.9 23.1 25.1 29.5 照射后 22.5 23.2 20.6 23.4 25.4 20.4 20.6 21.9 22.6 23.5 24.3 R程序如下：

30、>x<-c(25.7,24.4,21.1,25.2,26.4,23.8,21.5,22.9,23.1,25.1,29.5) > y<-c(22.5,23.2,20.6,23.4,25.4,20.4,20.6,21.9,22.6,23.5,24.3) > t.test(x,y,alternative="two.sided",paired=TRUE, conf.level=0.95) 結(jié)果如下： Paired t-test data: x and y t = 4.1334, df = 10, p-value = 0.002033 alternati

31、ve hypothesis: true difference in means is not equal to 0 95 percent confidence interval:（注：95%置信度樣本差值的可信區(qū)間） 0.850649 2.840260 sample estimates: mean of the differences 1.845455 結(jié)果大家可以自己分析。 例2：某豬場從10窩大白豬的仔豬中，每窩抽出性別相同、體重接近的仔豬2頭，將每窩兩頭仔豬隨機地分配到兩個飼料組，進行飼料對比試驗，試驗時間30天，增重結(jié)果見下表。試檢驗兩種飼料喂飼的仔豬平均增重差異是否顯著？ 窩號 1

32、2 3 4 5 6 7 8 9 10 飼料 10.0 11.2 12.1 10.5 11.1 9.8 10.8 12.5 12.0 9.9 飼料 9.5 10.5 11.8 9.5 12.0 8.8 9.7 11.2 11.0 9.0 R程序如下： > x<-c(10,11.2,12.1,10.5,11.1,9.8,10.8,12.5,12.0,9.9) >y<-c(9.5,10.5,11.8,9.5,12,8.8,9.7,11.2,11,9) > t.test(x,y,alternative="two.sided",paired=TRUE,c

33、onf.level=0.95) 結(jié)果如下： Paired t-test data: x and y t = 3.4553, df = 9, p-value = 0.007214 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.2382603 1.1417397 sample estimates: mean of the differences 0.69 三、兩個樣本的T檢驗 例1：分別測定了10只大耳白家兔、11只青紫藍家兔在停食18小時后正

34、常血糖值如下，問該兩個品種家兔的正常血糖值是否有顯著差異？（單位：kg） 大耳白 57 120 101 137 119 117 104 73 53 68 青紫藍 89 36 82 50 39 32 57 82 96 31 88 R程序如下： > x<-c(57,120,101,137,119,117,104,73,53,68) > y<-c(89,36,82,50,39,32,57,82,96,31,88) 首先進行方差齊次性分析。 > var.test(x,y) F test to compare two variances data: x and y F =

35、1.3454, num df = 9, denom df = 10, p-value = 0.6479 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.3560266 5.3330292 sample estimates: ratio of variances 1.345411 從結(jié)果中可以看出，方差可以認為是同質(zhì)的。 > t.test(x,y,alternative="two.sided",paired=FALSE,

36、var.equal=TRUE,conf.level=0.95) Two Sample t-test data: x and y t = 2.7179, df = 19, p-value = 0.01365 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 7.564314 58.235686 sample estimates: mean of x mean of y 94.9 62.0 例2：有人曾對公雛雞作了性激素效應試驗。將22只公雛雞完全隨機

37、地分為兩組，每組11只。一組接受性激素A（睪丸激素）處理；另一組接受激素C（雄甾烯醇酮）處理。在第15天取它們的雞冠個別稱重，所得數(shù)據(jù)如下： 激素 雞冠重量（mg）A 57 120 101 137 119 117 104 73 53 68 118 C 89 30 82 50 39 22 57 32 96 31 88 R程序如下： > x<-c(57,120,101,137,119,117,104,73,53,68,118) > y<-c(89,30,82,50,39,22,57,32,96,31,88) > var.test(x,y) F test to comp

38、are two variances data: x and y F = 1.0934, num df = 10, denom df = 10, p-value = 0.8904 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.2941901 4.0641018 sample estimates: ratio of variances 1.093443 >t.test(x,y,alternative="two.sided"

39、;,paired=FALSE,var.equal=TRUE,conf.level=0.95) Two Sample t-test data: x and y t = 3.3764, df = 20, p-value = 0.003 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 15.66997 66.33003 sample estimates: mean of x mean of y 97 56 第四節(jié) 方差分析一、 完全隨機資料的方差分析

40、方法1：采用rep函數(shù)和factor函數(shù) 例： A 組：12,15,21,31,21,15 B 組：21,22,15,15,21 C 組：24,15,12,15,32,15,12 程序如下： > c(12,15,21,31,21,15,21,22,15,15,21,24,15,12,15,32,15,12)->x > factor(c(rep(1,6),rep(2,5),rep(3,7)->group > anova(lm(xgroup) 注：1.rep(1,6)函數(shù)表示：將1重復6次。 2.factor()函數(shù)則是表明group是一個分組變量 3.lm(xgro

41、up)表示建立x為因變量，group為自變量的線性模型 4.anova(lm(xgroup)表示將線性模型中的方差表輸出。 結(jié)果如下： Analysis of Variance Table Response: x Df Sum Sq Mean Sq F value Pr(>F) group 2 5.95 2.98 0.0729 0.93 Residuals 15 612.49 40.83 方法2 例：某水產(chǎn)研究所為了比較四種不同配合飼料對魚的飼喂效果，選取了條件基本相同的魚20尾，隨機分成四組，投喂不同飼料，經(jīng)一個月試驗以后，各組魚的增重結(jié)果列于下表。 飼料 魚的增重（xij） 合計.i

42、x 平均.ix A1 31.9 27.9 31.8 28.4 35.9 155.9 31.18 A2 24.8 25.7 26.8 27.9 26.2 131.4 26.28 A3 22.1 23.6 27.3 24.9 25.8 123.7 24.74 A4 27.0 30.8 29.0 24.5 28.5 139.8 27.96合計 .x=550.8 R語言程序如下：建立一個文本文件保存在硬盤上,內(nèi)容如下: weight group31.9 1 27.9 1 31.8 1 28.4 135.9 1 24.8 2 25.7 2 26.8 2 27.9 2 26.2 2 22.1 3 23.6

43、 3 27.3 3 24.9 3 25.8 3 27.0 4 30.8 4 29.0 4 24.5 428.5 4 （用group來代表分組） > read.table("c:/r/6-1.txt",header=T)->test > edit(test)->test1 （這一步將group列的類型修改為字符,然后保存在test1向量中）. > aov(weightgroup,data=test1)->z > summary(z) 結(jié)果如下： Df Sum Sq Mean Sq F value Pr(>F) group 3 1

44、14.268 38.089 7.1362 0.002942 * Residuals 16 85.400 5.337 - Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1 二、兩因素方差分析 方法1：采用rep函數(shù)和factor函數(shù) 例：為了比較4種飼料(A)和豬的3個品種(B)，從每個品種隨機抽取4頭豬(共12頭)分別喂以4種不同飼料。隨機配置，分欄飼養(yǎng)、位置隨機排列。從60日齡起到90日齡的時期內(nèi)分別測出每頭豬的日增重(g),數(shù)據(jù)如下，試檢驗飼料

45、及品種間的差異顯著性。(FA=11.13，F(xiàn)B=13.21，MSe=202.0833)。 A1 A2 A3 A4 B1 505 545 590 530 B2 490 515 535 505 B3 445 515 510 495 R程序如下：x<-scan() 1: 505 545 590 530 490 515 535 505 445 515 510 495 13: Read 12 items > g1<-factor(c(rep(1:4,3) > g2<-factor(c(rep(1,4),rep(2,4),rep(3,4) > anova(lm(xg1+

46、g2) 結(jié)果如下： Analysis of Variance Table Response: x Df Sum Sq Mean Sq F value Pr(>F)g1 3 6750.0 2250.0 11.134 0.007269 * g2 2 5337.5 2668.8 13.206 0.006343 * Residuals 6 1212.5 202.1 - Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1 方法2： 例：為研究雌激素對子宮

47、發(fā)育的影響，現(xiàn)有4窩不同品系未成年的大白鼠，每窩3只，隨機分別注射不同劑量的雌激素，然后在相同條件下試驗，并稱得它們的子宮重量，見表6-21，試作方差分析。 6-21 各品系大白鼠不同劑量雌激素的子宮重量(g) 雌激素注射劑量(mg/100g)(B) 品系(A) B1(0.2) B2(0.4) B3(0.8) 合計xi. 平均.ix A1 106 116 145 367 122.3 A2 42 68 115 225 75.0 A3 70 111 133 314 104.7 A4 42 63 87 192 64.0 合計x.j 260 358 480 1098 平均jx. 65.0 89.5 1

48、20.0 R程序如下： 在硬盤上建立一個文本文件，內(nèi)容如下： weight ga gb 106 1 1 116 1 2 145 1 3 42 2 1 68 2 2 115 2 3 70 3 1 111 3 2 133 3 3 42 4 1 63 4 287 4 3 保存。 > read.table("c:/r/6.5.txt",header=T)->test > edit(test)->test1 > anova(lm(weightga+gb,data=test1)->z > z Analysis of Variance Table

49、 Response: weight Df Sum Sq Mean Sq F value Pr(>F) ga 3 6457.7 2152.6 23.771 0.0009923 * gb 2 6074.0 3037.0 33.537 0.0005535 * Residuals 6 543.3 90.6 - 三、有重復觀測值的兩因素方差分析和交互作用分析 沒有重復觀測值的數(shù)據(jù)，不能進行交互作用的方法分析。這是因為： (1)在這種情況下，(6-31)式中SSe,dfe實際上是A、B兩因素交互作用平方和與自由度，所算得的MSe是交互作用均方，主要反映由交互作用引起的變異。(2)這時若仍按上面所采用

50、的方法進行方差分析，由于誤差均方值大(包含交互作用在內(nèi))，有可能掩蓋試驗因素的顯著性，從而增大犯型錯誤的概率。(3)因為每個水平組合只有一個觀測值，所以無法估計真正的試驗誤差，因而不可能對因素的交互作用進行研究。 因此，進行兩因素或多因素試驗時，一般應設置重復，以便正確估計試驗誤差，深入研究因素間的交互作用。 例：為了研究飼料中鈣磷含量對幼豬生長發(fā)育的影響，將鈣(A)、磷(B)在飼料中的含量各分4個水平進行交叉分組試驗。先用品種、性別、日齡相同，初始體重基本一致的幼豬48頭，隨機分成16組，每組3頭，用能量、蛋白質(zhì)含量相同的飼料在不同鈣磷用量搭配下各喂一組豬，經(jīng)兩月試驗，幼豬增重結(jié)果(kg)列

51、于表6-29，試分析鈣磷對幼豬生長發(fā)育的影響。 表6-29 不同鈣磷用量(%)的試驗豬增重結(jié)果(kg) B1(0.8) B2(0.6) B3(0.4) B4(0.2) Ai合計xi. Ai平均.ix x1j1 22.0 30.0 32.4 30.5 A1(1.0) 26.5 27.5 26.5 27.0 24.4 26.0 27.0 25.1 324.9 27.1 x1jl 72.9 83.5 85.9 82.6 1jx 24.3 27.8 28.6 27.5 x2j. 23.5 33.2 38.0 26.5A2(0.8) 25.8 28.5 35.5 24.0 x2jl 27.0 30.1

52、33.0 25.0 350.1 29.2 76.3 91.8 106.5 75.5 x2jx 25.4 30.6 35.5 25.2 x3jl 30.5 36.5 28.0 20.5 26.8 34.0 30.5 22.5 25.5 33.5 24.6 19.5 332.4 27.7 x3j. 82.8 104.0 83.1 62.5 A3(0.6) x3jx 27.6 34.7 27.7 20.8 34.5 29.0 27.5 18.531.4 27.5 26.3 20.0 x4jl 29.3 28.0 28.5 19.0 319.5 26.6 x4j. 95.2 84.5 82.3 57.

53、5 A4(0.4) . 4jx 31.7 28.2 27.4 19.2 Bj合計x.j. 327.2 363.8 357.8 278.1 1326.9 Bj平均.jx 27.3 30.3 29.8 23.2 27.6 R程序如下： 建立一個文本文件，內(nèi)容和格式如下： weight ga gb 22 1 1 26.5 1 124.4 1 1 30 1 2 27.5 1 2 26 1 2 32.4 1 3 26.5 1 3 27 1 3 30.5 1 4 27 1 4 27.5 4 2 28 4 2 27.5 4 3 26.3 4 3 28.5 4 3 18.5 4 4 20 4 4 19 4 4

54、 將該文件保存在硬盤上。 > read.table("c:/r/6-6.txt",header=T)->test > edit(test)->test1 > anova(lm(weightga+gb+ga:gb,data=test1)->su > su 結(jié)果如下： Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) ga 3 44.51 14.84 3.2207 0.03558 * gb 3 383.74 127.91 27.

55、7667 4.919e-09 * ga:gb 9 406.66 45.18 9.8085 5.115e-07 * Residuals 32 147.41 4.61 - Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1 四、多重比較 多重比較需要安裝兩個package，一個是mvtnorm，一個是multcomp。 > library(mvtnorm) > library(multcomp) 1 單因素方差分析多重比較 例：某水產(chǎn)研究所為了比較四種不同配合飼料對魚的飼喂效果，選取了條件基本相同的魚20尾，隨機分成四組，投喂不