R語言及在農業試驗數據分析中的基本應用（二）

王磊 閔佳鑫 申紅芳 鄂志國

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

在R 語言[1]及在農業試驗數據分析基本應用的兩篇文章的第一篇[2]，對R 語言作了初步介紹以及非常基本的簡單應用后，本文采用蓋鈞鎰主編的教材《試驗統計方法》[3]中的示例數據集，重點介紹利用R 語言在農業試驗數據分析中的一些較為常用的方法，主要有t檢驗、方差分析以及均值的多重比較、卡平方檢驗等，最后提出R 語言使用的幾點建議。

1 兩個樣本平均數的比較

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

1.1 成組數據的平均數比較

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

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

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

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

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

1.2 成對數據的平均數比較

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

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

表1 A、B 兩處理病毒在番茄上產生的病痕數

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

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

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

2 方差分析

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

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

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

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

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

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

表3 例3 的4 種藥劑的試驗數據的方差分析結果

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

3 均值多重比較

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

3.1 LSD 法

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

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

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

或者，直接寫為：

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

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

3.2 Duncan 法

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

3.3 Tukey 法

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

3.4 SNK 法

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

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

4 卡平方檢驗

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

4.1 適合性卡方檢驗

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

表4 F2 代表型的觀察次數

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

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

4.2 列聯表獨立性卡方檢驗

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

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

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

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

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

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

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

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

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

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

5 總結和建議

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

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

在利用R分析數據時，有幾點建議：

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

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

圖1 例3數據用以讀入的csv 的格式

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

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

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

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

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

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

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

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

data：用于設置分析的數據集。

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

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

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