R語言及在農(nóng)業(yè)試驗數(shù)據(jù)分析中的基本應用（二）

王磊 閔佳鑫 申紅芳 鄂志國

（中國水稻研究所/水稻生物學國家重點實驗室，杭州 310006；第一作者：wanglei05@caas.cn；*通訊作者：ezhiguo@caas.cn）

在R 語言[1]及在農(nóng)業(yè)試驗數(shù)據(jù)分析基本應用的兩篇文章的第一篇[2]，對R 語言作了初步介紹以及非?；镜暮唵螒煤?，本文采用蓋鈞鎰主編的教材《試驗統(tǒng)計方法》[3]中的示例數(shù)據(jù)集，重點介紹利用R 語言在農(nóng)業(yè)試驗數(shù)據(jù)分析中的一些較為常用的方法，主要有t檢驗、方差分析以及均值的多重比較、卡平方檢驗等，最后提出R 語言使用的幾點建議。

1 兩個樣本平均數(shù)的比較

假設有兩組樣本，我們希望通過t 檢驗比較這兩組樣本平均數(shù)，以檢驗兩組樣本所屬的總體均值有無差異。因為兩組樣本數(shù)據(jù)對應的試驗設計和取樣等方面的不同而分為成組數(shù)據(jù)的平均數(shù)比較和成對數(shù)據(jù)平均值比較兩種情形。t 檢驗用到的函數(shù)是基礎程序包stats 中的函數(shù)t.test()。我們以《試驗統(tǒng)計方法》中的兩個示例作為例子，介紹利用t.test()作出成組數(shù)據(jù)和成對數(shù)據(jù)的t 檢驗。

1.1 成組數(shù)據(jù)的平均數(shù)比較

我們以《試驗統(tǒng)計方法》中的例5.3 數(shù)據(jù)集[3]作為例子，介紹利用函數(shù)t.test()進行成組數(shù)據(jù)的平均數(shù)比較。

例1 兩種密度產(chǎn)量的差異。調(diào)查某農(nóng)場栽插密度為30萬苗/667 m2和35萬苗/667 m2的稻田各5 塊，得平均產(chǎn)量（kg/667 m2）：400、420、435、460、425和450、440、445、445、420，試檢驗兩種密度單位面積（667 m2）產(chǎn)量的差異顯著性。

根據(jù)題意，我們希望檢驗兩種密度單位面積（667 m2）產(chǎn)量對應的總體均值的差異顯著性。利用函數(shù)c()分別輸入密度為30萬苗和35萬苗5 塊稻田的單位面積（667 m2）產(chǎn)量：

從計算結(jié)果可知，兩組樣本數(shù)據(jù)的平均值分別是428 kg/667 m2和440/667 m2，t值=－1.0776，自由度（df）=6.1086，P值（p－value）=0.3219，大于0.05。另外，我們注意到計算結(jié)果的標題是Welch Two Sample t－test，即兩組樣本的Welch t 檢驗，這是兩組樣本對應的總體方差不要求相等的t 檢驗。函數(shù)t.test()對兩個總體方差的默認設置是不等：var.equal=FALSE，對應的t 檢驗即為Welch t 檢驗。如果假設兩樣本的總體方差相同，對應的函數(shù)參數(shù)var.equal 的設置為var.equal=TRUE，這是常規(guī)的t 檢驗，這時函數(shù)t.test()的t 檢驗為：

我們注意到常規(guī)的t 檢驗的P值=0.3126，略小于Welch t 檢驗得到的P值。由于Welch t 檢驗不要求兩組樣本對應的總體方差相等，所以Welch t 檢驗的結(jié)果更為可靠，在實際的t 檢驗中，推薦使用Welch t 檢驗。對于這一例子的數(shù)據(jù)，30萬苗/667 m2和35萬苗/667 m2的5 塊稻田單位面積（667 m2）平均產(chǎn)量相差12 kg，差異不大，基于常規(guī)的t 檢驗和Welch t 檢驗，統(tǒng)計上都不顯著（P＞0.05）。

1.2 成對數(shù)據(jù)的平均數(shù)比較

我們以《試驗統(tǒng)計方法》中例5.6 的數(shù)據(jù)集[3]作為例子介紹成對樣本數(shù)據(jù)t 檢驗，這時，兩組數(shù)據(jù)的總體方差是否相等已經(jīng)無關緊要了。

例2 選生長期、發(fā)育進度、植株大小和其他方面都比較一致的兩株番茄構(gòu)成一組，共得7 組，每組中一株接種A 處理病毒，另一株接種B 處理病毒，以研究A處理和B 處理方法的鈍化病毒效果。處理A和處理B病毒在番茄上產(chǎn)生的病痕數(shù)目見表1，試檢驗兩種處理方法的差異顯著性。

表1 A、B 兩處理病毒在番茄上產(chǎn)生的病痕數(shù)

成對樣本數(shù)據(jù)的檢驗可以采用兩種方法。第一種方法是類似于成組樣本數(shù)據(jù)的t 檢驗，但利用函數(shù)t.test()的設置兩組數(shù)據(jù)是否成對的參數(shù)paired，paired 默認設置是否，即paired=FALSE，因為是成對數(shù)據(jù)，所以paired 需要設置為真，即paired=TRUE。

還有一種方法，是先計算成對數(shù)據(jù)之差d，然后利用函數(shù)t.test()對計算得到的差數(shù)d 進行單樣本的t 檢驗。

顯然兩種計算方法計算結(jié)果相同。根據(jù)計算結(jié)果，我們可知，差數(shù)平均值（mean of the difference 或者mean of x）= －0.8285714，即處理A和B 對應的病毒在番茄上產(chǎn)生平均病痕數(shù)之差約為－8.3，t 檢驗的t值=－4.1499，自由度（df）= 6，對應的P值（p－value）=0.006012 ＜0.01，所以兩種病毒接種處理方法的效果差異極顯著。由此我們可以得出結(jié)論，A 處理病毒與B 處理病毒的鈍化效果產(chǎn)生的病痕數(shù)平均值差異是－8.3，差異頗大，而且統(tǒng)計上也極顯著。

2 方差分析

比較兩組數(shù)據(jù)均值用的是t 檢驗，比較兩組以上的均值需要采用方差分析方法。方差分析的計算可以采用基礎包stats 中的函數(shù)aov()或者lm()。下面以《試驗統(tǒng)計方法》中的一個單向分組數(shù)據(jù)示例（例6.1和例6.3）[3]介紹用函數(shù)aov()對單向分組數(shù)據(jù)的方差分析。其他試驗設計類型數(shù)據(jù)的方差分析大致類似，只是不同的試驗設計對應函數(shù)aov()中不同的模型公式。

例3 以A、B、C、D 4 種藥劑處理水稻種子，其中A 為對照，每個處理各得4個苗高觀察值（cm），其結(jié)果如表2 。試對試驗數(shù)據(jù)進行分析。

表2 水稻施用不同藥劑處理的苗高（cm）

首先利用函數(shù)c()按照藥劑從A 到D 分別依次錄入各自的4個苗高數(shù)據(jù)，保存為數(shù)據(jù)向量y，而相應的藥劑分類變量是用向量元素重復函數(shù)rep()生成的，保存為字符向量group，其中LETTERS 是R 的內(nèi)置大寫字母常量，A、B、C、D 只需用LETTERS[1:4]表達即可，each=4 表示4個字母分別重復4次，而字符#表示注釋行標記，該字符的右側(cè)的注釋內(nèi)容代碼運行時跳過忽略。

我們利用函數(shù)aov()進行方差分析，將方差分析結(jié)果保存為myaov，然后利用函數(shù)summary()給出方差分析表。aov()中的設置y ～group 是單向分組數(shù)據(jù)對應的方差分析模型公式。不同的試驗設計或者數(shù)據(jù)結(jié)構(gòu)（如雙向分組數(shù)據(jù)等）需要采用各自特定的方差分析模型公式。

基于計算得到的方差分析表，我們給出了標準的方差分析表（表3），表中也列出了計算結(jié)果的英文名稱。在論文寫作時，欄目表頭的英文術語有所不同，如Sum Sq和Mean Sq 一般寫為SS和MS。計算結(jié)果中的最后一行Signif. codes 表示顯著性標記的說明：如果顯著性F 檢驗的P值＜0.001，用***標記；P值＜0.01，用** 標記；P值＜0.05，用* 標記；P值＜0.1，用.標記；P值＞0.1，不標記。同時也注意到，在計算結(jié)果中并沒有給出方差分析表（表3）中最后一行“總計（Total）”中的“自由度”和“平方和”的數(shù)值，這兩個數(shù)值需要我們自己計算。

表3 例3 的4 種藥劑的試驗數(shù)據(jù)的方差分析結(jié)果

在計算得到的方差分析表中的第一列最右邊的數(shù)字是藥劑因子group 的F 檢驗對應的P值，等于5.10e－05，顯然小于0.001（右邊標注為***），4 種藥劑相互之間的差異極其顯著（P ＜0.001），說明試驗數(shù)據(jù)有非常強的證據(jù)表明4 種藥劑中至少有一對總體均值是不等的。接下來，我們需要探究4 種藥劑中兩兩之間的均值差異的顯著性。

3 均值多重比較

因為例3 方差分析中的藥劑主效極顯著，不同藥劑對水稻苗高有不同效應，那么試驗者自然會感興趣到底是其中那幾對藥劑之間有不同效應，這就需要我們對這4 種藥劑的總體均值之間的差異性進行統(tǒng)計檢驗，即開展藥劑均值之間的多重比較。常用的多重比較方法有最小顯著差數(shù)檢驗法（LSD 法）、Duncan 新復極差法（Duncan 法）、Tukey 固定極差檢驗法（Tukey 法）和Student－Newman－Keul 復極差檢驗法（SNK 法或NK法）。不少程序包都提供了這些多重比較的檢驗方法，這里我們采用程序包agricolae 的函數(shù)LSD.test()、duncan.test()、HSD.test()和SNK.test()進行相應的多重比較[4]，這幾個函數(shù)來自同一程序包，使用方法相同。

3.1 LSD 法

在例3 的方差分析中，我們已經(jīng)利用函數(shù)aov()對試驗數(shù)據(jù)作了方差分析，基于函數(shù)aov()的計算結(jié)果保存為myaov，我們將采用程序包agricolae 的函數(shù)LSD.test()對藥劑的藥效進行LSD 多重比較。首先需要裝載程序包agricolae，然后再調(diào)用函數(shù)LSD.test()，其中藥劑變量的名稱為group（兩側(cè)需添加雙引號），alpha=0.05表示顯著性水平α = 0.05，console=TRUE 表示在控制臺窗口顯示計算結(jié)果（默認為不顯示）。

輸出結(jié)果較多，我們注意到以字母表示的多重比較結(jié)果是最后一部分。為了節(jié)省篇幅，我們只給出字母表示的多重比較結(jié)果，為此，我們先將多重比較的計算結(jié)果保存為LSD5，然后利用函數(shù)names()查看輸出結(jié)果的組件：

所以，LSD5 共有5個組件部分，其中groups 是用字母表示的多重比較排列結(jié)果，我們利用美元$符號獲取該部分的輸出結(jié)果：

或者，直接寫為：

注意，此時用以設置是否在控制臺窗口顯示計算結(jié)果的參數(shù)console 不需設置為TRUE，默認設置即可，因而函數(shù)中不必出現(xiàn)該參數(shù)。

類似地，我們可以計算得出，1%顯著性水平下的多重比較結(jié)果：

3.2 Duncan 法

Duncan 新復極差法的多重比較可以采用程序包agricolae 中的函數(shù)duncan.test()。duncan.test()的用法與LSD.test()用法相同，得到5%和1%顯著性水平的Duncan 新復極差多重比較的結(jié)果。

3.3 Tukey 法

這是TUKEY 在1952年提出的一種多重比較方法，該方法以控制試驗錯誤率為目標，又叫固定極差的q 檢驗法[5]。我們采用agricolae 中的函數(shù)HSD.test()進行Tukey 法的多重比較。用法與LSD.test()和duncan.test()相同。

3.4 SNK 法

SNK 法也被稱為q 法[3]，它是Tukey 法的一個發(fā)展，相對于Tukey 法，顯得較為不保守（傾向于發(fā)現(xiàn)較多的差異）。我們采用程序包agricolae 的函數(shù)SNK.test()比較4 種藥劑均值，類似于前述的幾種多種比較方法，得到5%和1%顯著性水平的SNK 法多重比較的結(jié)果。

從計算可知，5%水平下4 種不同方法的多重比較的結(jié)果有所不同，但在1%水平下，結(jié)果相同，效果最好的藥劑D 極顯著地好于效果較差的藥劑A和C，藥劑A和C 沒有顯著差異，效果排第2 的藥劑B 只是與藥劑C 有極顯著差異。

4 卡平方檢驗

卡平方（χ2）檢驗可用于樣本間的方差同質(zhì)性比較、計數(shù)數(shù)據(jù)的適合性以及基于列聯(lián)表的兩個變量的獨立性等問題。我們利用基礎包stats 中的函數(shù)chisq.test()對《試驗統(tǒng)計方法》中的1個兩對等位基因遺傳試驗的數(shù)據(jù)集（例7.6）[3]進行適合性檢驗以及1個列聯(lián)表（例7.9）進行列變量和行變量的獨立性檢驗。

4.1 適合性卡方檢驗

例4 兩對等位基因遺傳試驗，如果基因為獨立分配，則F2代4 種表現(xiàn)型在理論上的比率為9∶3∶3∶1。有一水稻遺傳試驗，以稃尖有色非糯品種與稃尖無色糯性品種雜交。其F2代品種得表4 結(jié)果。試檢查實際結(jié)果是否符合為9∶3∶3∶1 的理論比率。

表4 F2 代表型的觀察次數(shù)

我們用函數(shù)chisq.test()作試驗數(shù)據(jù)與理論比率的適合性卡方檢驗，其中x和p 是函數(shù)參數(shù)，用于設置試驗觀察值頻數(shù)以及理論比率。

從結(jié)果可知卡方值（X－squared）= 92.706，P值（pvalue）＜2.2e－16，幾乎為0，從而試驗數(shù)據(jù)有極強的證據(jù)拒絕原假設，即有極強的證據(jù)表明該水稻稃尖和糯性性狀在F2的稃尖表型分類結(jié)果不符合9∶3∶3∶1 的理論比率，也就是說，該兩對等位基因并非獨立遺傳，而可能為連鎖遺傳。

4.2 列聯(lián)表獨立性卡方檢驗

例5 進行大豆等位酶Aph 的電泳分析，193 份野生大豆和223 份栽培大豆的等位基因型的次數(shù)列于表5 中，試分析大豆Aph等位酶的等位基因型頻率是否因物種而不同。

表5 野生大豆和栽培大豆Aph等位酶的等位基因型次數(shù)分布

對于試驗得到的列聯(lián)表數(shù)據(jù)，研究者感興趣檢驗大豆等位酶Aph 的等位基因型頻率是否因物種而不同，對應的統(tǒng)計檢驗是獨立性檢驗，相應的原假設H0和備選假設Ha分別是：

H0：大豆Aph等位酶的等位基因型頻率與物種無關；

Ha：大豆Aph等位酶的等位基因型頻率與物種有關。

我們?nèi)匀焕煤瘮?shù)chisq.test()進行列聯(lián)表的獨立性檢驗，但列聯(lián)表數(shù)據(jù)必須以矩陣的格式輸入，用到的函數(shù)是矩陣函數(shù)matrix()，其中用到2個參數(shù)：第一個參試nrow 用以設置矩陣行數(shù)，這里矩陣行數(shù)為2，所以nrow=2；另外一個參數(shù)是byrow 用以設置數(shù)據(jù)在矩陣的排列順序，默認值為FALSE，按照列的順序排列，即byrow=FALSE。

我們按照參數(shù)byrow 的默認設置byrow=FALSE和按照行的順序排列的設置byrow=TRUE 分別輸入數(shù)據(jù)生成矩陣mymat，體會兩種不同設置的差別（注意其中函數(shù)c()中的數(shù)據(jù)的不同排列順序）。

# 按照默認的列的順序排列

然后利用函數(shù)chisq.test()進行列聯(lián)表的獨立性檢驗：

從計算結(jié)果可知，卡方值（X－squared）= 154.04,自由度（df）= 2, 相應的P值（p－value）＜2.2e－16 ，所以數(shù)據(jù)有非常強的證據(jù)拒絕H0，即野生大豆和栽培大豆的Aph等位基因型頻率有顯著差別。

5 總結(jié)和建議

我們介紹了如何利用R 語言對試驗數(shù)據(jù)進行t 檢驗、方差分析、均值的多重比較以及卡平方檢驗等常用的統(tǒng)計分析方法，在分析計算中，只是需要掌握加載已經(jīng)下載安裝的所需程序包，然后調(diào)用程序包中的函數(shù)，設置好函數(shù)相關的參數(shù)即可，并不需要關心具體的統(tǒng)計分析是如何做的。如果這些也算作編程的話，那么這樣的編程應該是比較容易掌握的。

不同類型和不同特點的數(shù)據(jù)、以及不同的研究目的需要不同的分析方法，而R 提供了不同研究領域上萬個多種多樣的程序包[6]，當然，其中有一些是一般性的統(tǒng)計分析程序包，如R 自帶的基礎程序包stats、混合線性模型分析程序包lme4等，不過，不同的研究領域都有對應的程序包可以選擇使用，如本文介紹的程序包agricolae 就是由秘魯科學家Felipe de Mendiburu 在農(nóng)業(yè)研究機構(gòu)工作時開發(fā)的專門用于農(nóng)業(yè)科研數(shù)據(jù)分析的程序包[4]。所以，R 語言就在一般的商業(yè)化軟件與專業(yè)程序員之間建起了一個豐富的折中選擇。這是R語言的特點，不僅體現(xiàn)了R 語言開發(fā)者研發(fā)的初衷，也是R 語言在世界范圍如此受歡迎的重要原因[7]。

在利用R分析數(shù)據(jù)時，有幾點建議：

1）數(shù)據(jù)的輸入或者讀入是分析的起點，也是關鍵點，其中尤其要注意R 函數(shù)所要求的數(shù)據(jù)格式大多是長形（long－format）。如對例3 的數(shù)據(jù)方差分析，不能按照表2 的格式輸入數(shù)據(jù)，第1 列是藥劑類型，第2 列到第5 列是苗高的4次重復觀察值，而是應該以2 列的格式（長形）輸入數(shù)據(jù)，其中第一列是藥劑類型，第二列是對應的苗高觀察值數(shù)據(jù)：

如果是從文本文件或者Excel 工作表讀入，那么數(shù)據(jù)的格式應該類似，如從Excel 電子表保存為csv 格式的數(shù)據(jù)形式見圖1。

圖1 例3數(shù)據(jù)用以讀入的csv 的格式

2）本文在利用R 函數(shù)的數(shù)據(jù)分析中，都沒有用到函數(shù)中用于設置數(shù)據(jù)分析數(shù)據(jù)集的參數(shù)data，這是因為我們在R 控制臺窗口直接輸入或者生成數(shù)據(jù)變量，R 函數(shù)可以直接調(diào)用，如果數(shù)據(jù)集是讀入的，那么在相應的R 函數(shù)中需要利用參數(shù)選項data 設置分析數(shù)據(jù)集。例如，我們在例2 中讀入圖1 中所示的csv 格式的數(shù)據(jù)（文件名為gjy61.csv，保存在當前工作文件夾中），那么讀入數(shù)據(jù)集以及數(shù)據(jù)的方差分析為：

另外還有幾種比較常用的關聯(lián)變量的方法，尤其是有些R 函數(shù)不提供設置數(shù)據(jù)集的參數(shù)選項時會用到，有利用函數(shù)attach()或者with()關聯(lián)數(shù)據(jù)集中變量的方法，還有更直接利用美元符號的方法。

3）本文分析示例數(shù)據(jù)用到的幾個不同的R 函數(shù)，其中函數(shù)的參數(shù)設置都比較少。R 函數(shù)一般都提供了較多的參數(shù)選項，從而我們在具體的數(shù)據(jù)分析中可以按照分析要求選取合適的參數(shù)進行設置。具體的函數(shù)參數(shù)選項可以通過R 系統(tǒng)的幫助功能或者直接在網(wǎng)上查詢。如函數(shù)t.test()使用方法為：

其中的參數(shù)paired=FALSE和var.equal=FALSE 在例1和例2 的數(shù)據(jù)分析時已經(jīng)用到了，更多的其他各項參數(shù)設置解釋如下：

x和y= NULL：x和y 都為數(shù)據(jù)向量，其中y =NULL 表示如果是單樣本t 檢驗，y可以忽略；

alternative：用于設置是雙尾假設（"two.sided"）還是單尾假設（"less" 或"greater"），如在例1 中，試驗者在試驗開始前已經(jīng)比較肯定密植能夠提高產(chǎn)量，試驗的目的是希望對此檢驗，那么alternative 的設置為alternative="less"，或者簡寫為alt="less"。注意默認是alt="two.sided"，即如果t 檢驗是雙尾檢驗，那么這一選項可以忽略；

mu=0：表示檢驗總體均值是否為0(單樣本) 或者兩個樣本總體均值是否相等，如在例1 中試驗者希望檢驗密植能否增產(chǎn)5 kg 以上，那么mu=－5（負數(shù)是因為x 是較低產(chǎn)量30萬苗/667 m2的稻田產(chǎn)量）；

conf.level=0.95：表示總體均值（單樣本）或者兩樣本均值差數(shù)的置信區(qū)間的置信度設置，默認置信度為95%；

data：用于設置分析的數(shù)據(jù)集。

函數(shù)t.test()的參數(shù)選項設置，其中的許多參數(shù)的默認設置我們很多時候都是接受的，如alternative="two.sided", mu=0, paired=FALSE, var.equal=FALSE, conf.level=0.95，那么函數(shù)t.test() 中都可以忽略，正如我們在例1 的數(shù)據(jù)分析中所做的。

4）網(wǎng)上有豐富的R 語言的資源，這是R 語言流行的另外一個重要原因[7]。例如，想用R 分析試驗數(shù)據(jù)，但覺得不知如何著手時，一種比較快捷有效的辦法是通過網(wǎng)上查詢類似試驗數(shù)據(jù)的例子，將例子中的R 代碼拷貝到R 軟件中程序腳本窗口，并將例子中的數(shù)據(jù)替換為自己的試驗數(shù)據(jù)，然后運行代碼進行分析。不過在具體的分析過程中，我們至少需要清楚采用哪種統(tǒng)計分析方法分析自己的試驗數(shù)據(jù)，并能正確解讀分析結(jié)果。例如，在分析例1 的試驗數(shù)據(jù)時，比較兩組樣本的均值所采用的方法是t 檢驗，那么我們可以在百度按照兩個關鍵詞“R 語言t 檢驗”查詢相關內(nèi)容和例子。

通過本文的介紹，期望有更多的農(nóng)業(yè)科技工作者下載開始使用這一具有強大的統(tǒng)計計算功能、便捷的數(shù)據(jù)可視化系統(tǒng)以及免費開源的R 語言。