R 语言及在农业试验数据分析中的基本应用(一)
2021-12-02王磊闵佳鑫鄂志国
王磊 闵佳鑫 鄂志国
(中国水稻研究所/水稻生物学国家重点实验室, 杭州310006;第一作者:wanglei05@caas.cn;*通讯作者:ezhiguo@caas.cn)
R 语言已经是工业界和学术界的主流和标准计算语言之一,可能已是科学研究中数据分析最流行的计算语言[1-4],更有学者认为R 语言是“免费、易学、最热门、最公开透明、最全面、最新最快的统计软件”[3]。但在中国知网(https://www.cnki.net/),选择年份从 2000 年至2020 年,按照关键词“R 语言”或者“R 软件”对我国几种主要农学类学术期刊进行全文检索,检索到《中国农业科学》100 篇,《作物学报》39 篇,而《中国水稻科学》只有6 篇,《杂交水稻》3 篇(其中1 篇是综述)。检索结果中,不排除只是在综述或者讨论中提到名称,但也有可能个别论文所作的试验数据分析采用了R 语言但没有注明而没有检索到。不过,总的来说,检索到的文章不少是组学研究和分子生物学研究领域的工作,农业试验数据分析,特别是水稻科研和基础的大田试验数据分析很少采用R 语言,其中的原因可能是:(1)农业科研工作者的软件使用习惯,譬如一直习惯采用SAS或者SPSS 等商用软件;(2)对R 语言了解不多或者很少接触;(3)也有可能是 R 语言的使用环境(R 软件)不同于菜单式的软件而需要编写代码的特点,妨碍了R语言的普及使用。
我们以两篇文章的篇幅,介绍R 语言的基本知识点以及在农业试验数据上的基本应用,希望有助于农业科研工作者开始尝试使用R 语言分析试验数据。本文先介绍R 的下载安装、起步使用以及R 的界面以及各项功能,接着介绍程序包的安装和调用,然后是如何直接输入试验数据或者从文本或者Excel 工作表读入,最后通过简单的示例介绍利用R 语言作一些简单的统计分析。
1 R 语言基本介绍
R 语言官方网站(https://www.r-project.org/)的第一句话是这样介绍R 语言的:“R is a free software environment for statistical computing and graphics”[5],即“R是一个用于统计计算和绘图的自由软件环境”。这句话表明了R 的三大特色:强大的统计计算功能、便捷的数据可视化系统以及免费开源。R 软件是R 语言的实现环境,借助于程序包,提供了一套完整的数据处理、统计分析和作图等功能,譬如:各种类型的作图、常用的统计检验、方差分析、线性和非线性建模、各种多元统计分析以及时间序列分析等等,使用者只需根据数据分析目的,清楚分析所需的图形类型、统计计算或者统计模型,调用相应的程序包以及其中的函数,便可高效进行数据的统计分析。
需要说明,R 软件是R 语言的实现环境,很多时候,我们不加以区分R 软件或者R 语言,有时直接简称为R。
1.1 下载与安装
R 软件的下载和安装可以按照下列步骤执行:
1)在 R 的官方网站(http://www.r-project.org),点击蓝色字体显示的超链接“download R”,给出的是“CRAN 镜像网站”(CRAN Mirrors)页面,其中按照国别给出了各自国家的镜像站。在给出的镜像网站清单中选择一个镜像站,我们一般是选择国内(China)的镜像站,譬如清华大学(TUNA Team, Tsinghua University)对应的镜像站域名:https://mirrors.tuna.tsinghua.edu.cn/CRAN/,点击进入该镜像站后,有适合3 个不同操作系统的下载选择:“download R for Linux”、“download R for(Mac) OS X”、“download R for Windows”,如果电脑是微软Windows 操作系统,点击“download R for Windows”,在出现的下载页面,点击“Base”,在随后的页面上点击“download R 4.0.5 for Windows”下载 R,其中的R 4.0.5 表明R 的版本号是4.0.5,这也是目前最新的版本号(2021 年 5 月 10 日)。
2)下载后,双击下载的安装文件R-4.0.5-win.exe开始安装。在安装过程中需要选择安装的语言,默认为中文。 这里的中文选择只是表明安装过程是中文界面。
3)如果按照默认设置,软件将安装在文件夹c:Program FilesRR-4.0.5,当然也可以选择不同的位置安装。
4)安装结束后,在桌面会生成一个R 启动命令的快捷键,真正的R 启动命令是位于目录c:Program FilesRR-4.0.5bin 中的命令文件Rgui.exe。
1.2 开始使用
只需双击桌面上的R 快捷键(或者点击开始菜单上的R 按钮)便可打开R 软件窗口(图1)。
R 软件窗口(或者直接称为“R 窗口”)界面比较简洁,由三部分组成:菜单栏、快捷键按钮以及称为R 控制台(R Console)的交互式命令对话窗口。R 控制台是R 窗口的主要界面,如果编写脚本或者作图,还会用到编辑器窗口和图形窗口。
R 窗口的菜单栏中具体栏目有文件、编辑、查看、其他、程序包、窗口、帮助;而菜单栏下方列出了8 个快捷键按钮,对应常用的一些命令功能,从第一个的“打开程序脚本”按钮到最后一个“打印”按钮。如果我们将鼠标的光标对准按钮,光标下方就会显示按钮的中文名称。
R 控制台是一个交互式命令对话窗口,命令提示符是大于符号“>”,在该提示符右侧闪动的光标位置表示在此处可以输入指令,回车即得结果。比如,简单的加法2+3,在输入命令后,按回车键执行:
在命令2+3 下方给出了计算结果5,类似地,各种命令可在R 命令行中输入,回车既得结果。
结果5 左侧有标记[1],该标记表示计算输出结果的第一行(本意是这一行的第一个字符是输出结果的第一个字符),不过,如果输出结果多于一行,那么第二行中括号的数字表示该行的第一个字符的序号,如第一行有5 个数字,那么第二行中括号中的数字为6,即,表示第二行的第一个数字为输出结果的第6 个数字;依次类推。如利用英文冒号:生成1 到30 的自然数:
从生成的自然数序列很容易明白最左侧的中括号中的数字所代表的含义。不过,结果显示的行的宽度与当时控制台窗口的宽度有关。
在R 命令行中交互式输入命令时,我们可以使用键盘中的上下方向键查看已输入命令的历史记录。根据需要,我们就可以选择一个之前输入过的命令,按回车键可重新执行该命令,当然也可以修改之前的命令并执行。
我们注意到启动R 软件后,在R 控制台窗口中“>”提示符上方的说明给出了R 的版本号(R version 4.0.5),并指出了R 是自由软件并提示可以用命令license()或者licence()获取版权相关的信息,可以用contributors()查看R 语言开发过程中的重要贡献者,用citation()可以获取R 软件的参考文献格式(包括latex 引用格式)。同时也提醒,特别是对于初学者,可以用命令demo()很粗略地了解R 的基本功能和能作出的一些基本图形。我们可以尝试输入命令:
输入该命令回车后给出了4 个基础程序包(base、graphics、grDevices 和stats)的演示示范函数。譬如程序包graphics 有6 个演示函数,如果想了解演示函数的演示内容,例如函数graphics(与程序包同名),则可输入命令,然后回车:
按照演示的提示,逐步按回车键,可以了解R 的一些图形功能以及作出的图形。
另外,也提示了可以用命令help()获取帮助信息,譬如希望查看t 检验函数t.test()的帮助信息,可以采用命令:
而help.start()给出了HTML 浏览器来查看帮助文件:
在显示的HTML 页面上给出了非常丰富的R 帮助信息的汇总页面。
如果我们希望退出R,可以单击软件窗口(RGui窗口)的右上角的叉号按钮,或者直接输入命令q()然后回车:
随后出现的对话框询问“是否保存工作空间映像?”,有 3 个选项:“是”、“否”和“取消”。如果选择“是”,2 个新文件:“.Rhistory”和“.Rdata”将会在工作目录下创建,所有先前输入的R 命令都会保存在“.Rhistory”,而所有的工作空间保存在“.RData”。注意这2 个文件都没有前缀。下一次启动R 时,通过这2 个文件,前一次保存的工作空间和输入的命令将自动恢复还原。
1.3 程序包
R 软件是R 语言的实现环境,具体来说,R 软件或许可以看作是由基于R 语言的许多常规和现代统计方法的运行环境,许多统计方法作为R 的基础运行环境的一部分是由称为程序包(package)的函数来提供的。软件R 下载安装时,已经自动安装了基础和常用的30个程序包,包括最基础的程序包base、stats 和 graphics,以及数据集包datasets。但还有许许多多丰富多彩、日新月异的数据分析程序包,用户可以很方便下载使用。目前程序包的数目已经过万[1],这是R 的使用方便、功能强大的最主要的一个原因。
一个程序包是由相关的函数、数据及文档组成。除了R 安装时经自动安装的30 个程序包,如果我们在数据分析时需要用到还没有安装的某个程序包中的函数,则我们需要先下载安装该程序包,然后加载该程序包,才能使用所需函数。当然,程序包只需下载安装1次即可,只是每次使用程序包中的函数之前,必须加载该程序包。在R 软件菜单中的“程序包”栏目提供了“安装程序包”和“加载程序包”选项,由此可以进行下载安装以及加载程序包。但更为方便的方法是利用命令进行下载安装和加载。
1.3.1 安装程序包
利用函数install.packages()下载安装程序包,函数括号内输入需要安装的包名,并在两边加引号(注意是在英文输入状态下的双引号)。以安装农业科研数据分析的一个程序包agricolae 为例,安装该程序包:
需要注意,R 对于英文字母的大小写敏感,所以要时刻注意区分字母的大小写。另外,也需注意区分中英文的标点符号,在R 语言中都是使用英文输入法状态下的标点符号。
1.3.2 加载程序包
一个程序包下载安装后,如果我们需要使用程序包中的函数,我们必须先加载该程序包,然后才能使用其中的函数。我们是利用函数library()来加载程序包的,函数括号内写上需要加载的程序包名称(可以不加双引号)。以加载程序包为例:
加载后,就可以使用该程序包的函数。如果希望了解程序包中的函数,同样可以利用命令library()来查看:
另外,每次R 启动时,自动加载了6 个基础程序包base、stats、graphics、grDevices、methods、utils 和 数 据 集包datasets,5 个基础程序包中的函数可以直接调用,不需重复加载这些程序包,如我们这里用到的2 个函数install.packages()和library()分别是基础包 utils 和 base的函数,不必考虑加载函数所属的程序包,可以拿来就用。
而基础包中的数据集包datasets 中含有不同研究领域不同类型的数据集104 个(版本4.2.0),不需加载,可以直接调用包中的数据集。包中的数据集中有不少是农学和生物学学科的数据集,在学习R 语言的具体统计分析方法时,可以通过数据集名称很方便地调用作为学习练习用的示例数据。如果希望查看包中的数据集,类似于查看程序包中的函数,利用命令library()查询:
在显示的数据集清单中,第4 个数据集是草本植物吸收的二氧化碳(Carbon Dioxide Uptake in Grass Plants)数据集CO2,如果我们希望了解数据集的内容,可以直接输入数据集名即可:
如果只是想查看数据集前几行,可以利用函数head()查看(n=3 表示只显示前3 行。这一参数选项可以简写为3):
给出的数据集中的第一行是各列的变量名称,随后给出了数据集中的观察值,最左侧的整数1、2、3 表示观察值对应行的序号。R 中变量名称可以是中文,但不建议这么去做,而观察值具体数值可以是中文。
如果进一步希望了解该数据集的背景,包括每个变量的定义等等,可以通过帮助命令查看:
帮助信息中,给出了数据集的描述(Description)、格式(Format)、采集数据集的试验内容(Details)、数据来源及相关文献(Source)以及基于该数据集的统计分析例子(Examples)。
1.4 工作目录
在开始使用R 分析数据之前,最好先确定工作目录。菜单中“文件”栏目中的“改变工作目录”可以查看当前的工作目录,并作出更改。我们也可利用命令getwd()查看当前工作目录,如需更改,可以利用函数命令setwd()进行更改。
这是当前工作目录,如要将工作目录改为其他目录,如D:MyR,那么通过命令setwd()作此改变:
注意在命令中的目录路径中是单斜杠,与Windows 中写法中的单斜杠方向正好相反。R 中另外一种目录路径的写法是双斜杠,方向与Windows 中写法中的单斜杠一致:
1.5 帮助
R 语言的介绍和使用不仅在网上很容易查询,也有许多中文书籍可以购买阅读。或许可以从阅读R 的官方网站中的中文帮助手册开始:例如,R 的综合档案网络平台CRAN(https://cran.r-project.org)有不少很好的帮助文档,其中在非英文部分(Non-English Documents)的中文名目(Chinese)下有3 份非常实用的中文文档(https://cran.r-project.org/other-docs.html):R reference card(R 参考卡片)、Frequently asked questions(R常见问题解答)以及Statistics and R Reading Notes(统计学与R 读书笔记)。
R 系统本身也提供便捷的帮助(R 窗口菜单中“帮助”栏目可以查询,只是帮助信息都是英文的)。想获取在线帮助,可调用函数。例如,我们希望简单了解程序包,先加载该程序包:
或者只是用问号“?”(注意是英文输入法状态下的问号):
如果希望了解函数包中的函数,可以利用函数library()查看:
如果希望了解某一个函数的使用方法,如计算均值的函数mean(),输入命令:
或者简单地用问号:
函数mean()是基础程序包base 的函数,R 启动时,基础程序包base 自动加载,所以可以直接寻求帮助。如果希望查看某个特定程序包的函数,必须确认该程序包已加载,不然需要先加载该程序包。如果我们希望查看agricolae 中的LSD 多重比较函数LSD.test(),需要先加载程序包再寻求帮助:
1.6 R 的编程及图形界面
在R 控制台交互窗口的命令行,输入命令,回车得到结果,简单直接。而且,通过键盘中的上下箭头键,可以调用已经输入过的命令和内容,不用重复输入。不过用了一段时间,慢慢熟悉了一些,可以学习创建程序脚本,执行的命令可以重复使用,提高计算效率。我们可以从R 窗口菜单中的“文件”栏目中的“新建程序脚本”生成R 编辑器窗口,然后在该窗口编写代码,利用R编辑器窗口中的“编辑”栏目中的选项运行执行编写好的代码。
如果觉得R 平台的综合功能以及使用界面不是很友好,则可以了解一些方便用户使用和编程的界面较为友好的图形界面软件(可参考“R 语言7 个免费的GUI 图 像 界 面 工 具 说 明 ”:https://blog.csdn.net/lhy55040817/article/details/8484883),并从中做出选择。目前使用最广的是Rstudio,建议在安装R 软件后,同时也安装 Rstudio。Rstudio 下载网站是:https://www.rstudio.com。
2 R 的数据输入和读入
学习任何一个数据统计分析软件,其中最关键的一个内容是掌握这一软件是如何输入或者读入数据的。R 的常用数据输入和读入大致可以分为以下几类:(1)直接输入;(2)基于文本文件的读入;(3)基于 Excel数据表的读入。每一类都有多种的方法,这里介绍的方法较为常用而且相对简单。
2.1 直接输入
如果试验数据较少,我们可以从控制台窗口直接输入数据。直接输入数据需要利用拼接函数c(),不同数值之间用英文逗号分开,如果是字符,字符两边需要用英文引号,同时,输入的数据(称为数据向量)一般会保存,便于随后的分析。在《试验统计方法》中的表3.7给出了某水稻杂种第二代植株米粒性状的分离情况,分离类别和频数分别是:红米非糯,96;红米糯稻,37;白米非糯,31,白米糯稻,15[6]。我们利用函数分别输入分离类别和频数,并保存为f2group 和f2freq:
如果想了解数据输入的变量内容,可以利用函数print()查看(或者直接输入变量名称。这里,两个命令在同一行,命令之间需要用英文分号隔开):
进一步,利用函数sum()可以很方便地计算这4 种类别频数的总数:
2.2 读入文本文件数据
假设我们需要读入以逗号分隔的csv 格式的文本数据文件datafile.csv(常常从Excel 数据文件将数据电子表保存为csv 格式文件),最常用的方法是用函数read.csv()读入(保存为 mydata):
其中选项header=TRUE 表示datafile.csv 的数据表有表头(类似于数据集CO2,每一列都带有变量名)。我们还可以用参数sep 来设置数据之间的分隔符号,默认为逗号,如果是空格符号,sep="";如果是制表符Tab键,sep=" "。假设是制表符作为分隔符号,那么利用函数读入datafile.csv 的设置为:
如果数据文件datafile.csv 不在工作目录中,那么我们需要给出该文件所在的目录路径,如文件在文件夹D:MyR,那么相应的读入命令修改为:
读入文本数据另外一个常用的函数是read.table(),用法也与read.csv()类似。
2.3 读入Excel 数据表
Excel 文件中的数据表可以利用openxlsx、readxl等程序包中的函数读入。我们以程序包openxlsx 中的函数 read.xlsx()为例简单介绍读取xlsx 格式的Excel 电子表数据。
首先先下载安装并且加载程序包openxlsx(其中#为注释符号,表明#右边的注释内容在运行代码时跳过忽略):
假如数据存储在Excel 文件MyExcelDataFile.xlsx中的第2 个工作表,那么可以利用包中的函数读入该数据电子表:
当然,我们也可通过工作表的名称来读入,比如工作表名称是sheet2,那么可以通过下面的代码读入该工作表:
注意,read.xlsx()的读入方式默认读入工作表表头,同时Excel 文件名以及其中的工作表可以是中文名称。
如果需要读入较老版本的Excel 文件(xls 格式),程序包readxl 提供了函数 read_xls(),该函数的用法与read.xlsx()类似。
3 小试一下
夸R 好,说R 行,但百闻不如一见,百见不如一试。我们在这一节简单用一下R,看R 是如何方便高效地进行数据汇总、作图,以及线性相关和线性回归计算。
3.1 数据汇总计算与作图
我们利用莫惠栋著作《农业试验统计》中的一个用于介绍线性回归分析的例子(例10.1)的数据集[7]来介绍R 的数据汇总和基本的作图功能,在下一小节介绍用R 函数进行线性相关和线性回归分析。
例1. 许多害虫的发生都和气象条件有一定关系。山东临沂测定10 年间(1964—1973 年)7 月下旬的温雨系数(雨量mm/平均温度C)和大豆第二代造桥虫发生量(每百株大豆上虫数)的关系如表1,试建立回归方程。
我们利用函数c()分别输入表1 中的温雨系数的数据以及对应的百株虫密度,并保存为TRcoef 和Dworm:
表1 温雨系数和大豆第二代造桥虫虫口密度
查看一下数据向量TRcoef 和Dworm 具体的数值:
对于数据向量TRcoef 和Dworm,可以利用R 函数很方便进行相关的汇总计算,如利用函数mean() 和sd() 计算平均值和标准差,或者利用函数summary () 对Dworm 进行简单的数据汇总(计算结果依次给出):
如果需要对数据向量的每个元素作同样的计算,如对向量Dworm 的每个元素作log 变换,不需要对每个元素分别做log 变换,只需要将函数log () 作用于Dworm 即可:
我们也可以很方便地作出常见的统计图形,如利用函数hist()和boxplot()作出数据向量Dworm 的直方图(图2)和箱线图(图 3):
如果需要作出百株虫密度Dworm 与温雨系数TR-coef 的散点图,我们可以采用函数plot ()(在参数设置中,用于x 轴变量的TRcoef 在前,用于y 轴变量的Dworm 在后):
从作出的图形(图4)可以帮助我们很清晰地看到百株虫密度Dworm 与温雨系数TRcoef 呈一定程度的负的线性相关。
进一步,如果我们希望完善图形,利用函数plot()参数设置很容易做到。比如,希望x 坐标轴和y 坐标轴的坐标说明不是变量名而是中文说明,那么可以利用函数的参数选项xlab 和ylab 进行设置。在下一小节图5 中可以看到添加坐标轴说明的效果):
plot () 的其他更多的参数选项和相关内容可以利用命令help(plot) 查看了解。
3.2 相关分析和线性回归
进一步,如果我们希望计算百株虫密度Dworm 与温雨系数TRcoef 之间的相关系数,并进行相关性检验,我们利用函数cor.test()进行计算和检验:
从计算结果的最后一行可知,百株虫密度Dworm与温雨系数TRcoef 的相关系数r = -0.9201,统计检验的 P 值(p-value)= 0.0001618,检验结果极显著,即数据有极强的证据表明两个变量之间的总体相关系数不等于 0。
从作出的2 个变量的散点图(图5)中,已经清楚百株虫密度与温雨系数呈一定程度的负的线性相关,进一步,我们希望建立基于温雨系数的百株虫密度的线性回归模型,利用函数 lm()和summary()可以得到拟合系数、统计检验以及决定系数R2:
从计算结果的回归系数(Coefficients)可知,简单线性回归的截距和斜率的拟合系数在Estimates 一栏给出,分别是179.212 和-14.110,其显著性检验的P 值在最后一栏(Pr(>|t|)给出,分别是 9.02e-07 和 0.000162,统计上都是极显著,而决定系数R2是在计算结果的最后一部分的第二行给出,Multiple R-squared:0.8466,即R2= 0.8466,也就是说,84.66%的大豆第二代造桥虫发生量的变异可以由7 月下旬的温雨系数解释。
最后,如果还希望在散点图(图4)上添加拟合的回归直线,这可以由函数abline()很方便地做到,其中myfit 是前面线性回归中保存的 lm (Dworm~TRcoef) 的计算结果。
4 总结和建议
相对于许多商用的统计分析软件,R 免费开源,优势明显。“便宜或者免费无好货”的说法,对于R 是不成立的。在R 开始应用的近 20 年,众人拾柴火焰高,开发的程序包无数,对于科研项目的数据,不管是来自田间试验还是实验室的实验,R 语言都有相应的程序包可以使用,满足分析的要求,而且很多时候同时有多个程序包和多种函数可以选择。
R 的一个明显缺陷是计算结果和帮助信息是用英文表述的,不过,一种较为方便的规避途径是通过百度、搜狗、必应等搜索引擎查询相关的内容,尤其是相关的例子,可以快速上手并帮助正确解读R 计算结果;另外一种较为系统的了解学习是通过相关的中文资料和书籍来熟悉常用函数的用法和结果的解读,除了文中提到的R 的综合档案网络平台CRAN 给出的3 份中文文档,参考文献中也列出了我们常常用来学习查阅的书籍,可供参考[3-4,8-10]。当然,在网上书店还有许许多多的相关书籍,可以选择适合自己的书来入门学习。
R 的使用如果能掌握一些简单的编程能力当然更好,不过大多数时候并不需要,只需要能正确调用及使用试验数据分析所需的程序包中的函数就可以。通过本文的简单介绍,希望能打消初学者对R 的畏难情绪,激发使用R 的兴趣,尽快入门,掌握R 使用的一些基本方法。像学习任何一个软件,活学活用、边用边学是最好的、也可能是最有效的方法,所以打开电脑,下载安装R,开始使用,用R 分析自己的试验数据吧。