R for Beginers_cn_2.0

R for Beginers_cn_2.0 - R for Beginners Chinese Edition 2.0...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: R for Beginners Chinese Edition 2.0 Emmanuel Paradis ´ Institut des Sciences de l’Evolution Universit´ Montpellier II e F-34095 Montpellier c´dex 05 e France E-mail: paradis@isem.univ-montp2.fr Co-translated by: XF Wang, YH Xie, JT Li and GH Ding 中文版说明 “R for beginners”是 一 本 公 认 的 经 典 手 册 , 非 常 适 合R的 初 学 者 。 英 文 原 版初著于2002年,而此稿是基于作者在2005年重新修订的第二版。 A Emmanuel Paradis博 士 为 本 稿 提 供 了 原 版 所 有L TEX源 文 件 。 翻 译 工 作 由 四 名 志 愿者 共 同 完 成 (Chap1–2: 王 学 枫 ;Chap3: 谢 益 辉 ;Chap4: 李 军 焘 ;Chap5–7: 丁 国 徽 ) 。 由 华 东 师 范 大 学 汤 银 才 老 师 负 责 本 文 档的 编 辑 校 订。北京大学李东风老师审阅了全稿并提出了大量宝贵意见。在此一并表示 衷心感谢! 编译仓促,差错难免,亟盼诸R友襄正。任何意见请通过pwxf@hotmail.com 联系我们。 译者 2006年4月 版权 c 2002, 2005, Emmanuel Paradis Permission is granted to make and distribute copies, either in part or in full and in any language, of this document on any support provided the above copyright notice is included in all copies. Permission is granted to translate this document, either in part or in full, in any language provided the above copyright notice is included. 目录 1 导言 1 2 基本原理与概念 2.1 基本原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 对象的产生,排列及删除 . . . . . . . . . . . . . . . . . . . . . 2.3 在线帮助 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 5 7 3 R的 数 据 操 作 3.1 对象 . . . . . . . . . . . . . . . . . . . . 3.2 在文件中读写数据 . . . . . . . . . . . . 3.3 存储数据 . . . . . . . . . . . . . . . . . 3.4 生成数据 . . . . . . . . . . . . . . . . . 3.4.1 规则序列 . . . . . . . . . . . . . 3.4.2 随机序列 . . . . . . . . . . . . . 3.5 使用对象 . . . . . . . . . . . . . . . . . 3.5.1 创建对象 . . . . . . . . . . . . . 3.5.2 对象的类型转换 . . . . . . . . . 3.5.3 运算符 . . . . . . . . . . . . . . . 3.5.4 访问一个对象的数值:下标系统 3.5.5 访问对象的名称 . . . . . . . . . 3.5.6 数据编辑器 . . . . . . . . . . . . 3.5.7 数学运算和一些简单的函数 . . . 3.5.8 矩阵计算 . . . . . . . . . . . . . 4 R绘 图 4.1 管理绘图 . . . . . . . . . 4.1.1 打开多个绘图设备 4.1.2 图形的分割 . . . . 4.2 绘图函数 . . . . . . . . . 4.3 低级绘图命令 . . . . . . . 4.4 绘图参数 . . . . . . . . . 4.5 一个实例 . . . . . . . . . 4.6 grid 和lattice 包 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 9 11 14 15 15 18 19 19 24 26 27 30 32 32 34 . . . . . . . . 37 37 37 38 41 42 44 45 49 5 R的 统 计 分 析 5.1 关于方差分析的一个简单例子 5.2 公式 . . . . . . . . . . . . . . 5.3 泛型函数 . . . . . . . . . . . 5.4 包 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 56 58 59 62 6 R编 程 实 践 65 6.1 循环和向量化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2 用R写程序 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.3 编写你自己的函数 . . . . . . . . . . . . . . . . . . . . . . . . . 68 7 R 相关的文献 72 1 导言 该手册是关于R的一个入门教材.由于主要针对初学者,我将重点放在了 对R的工作原理的解释上。R涉及广泛,因此对于初学者来讲,了解和掌握一 些基本概念及原理是很有必要的。在打下扎实的基础后,进行更深入的学习 将会变得轻松许多。本着深入浅出的宗旨,本手册将大量配合图表等形式, 尽可能使用通俗的语言,使读者容易理解而并不失细节。 R是一个有着统计分析功能及强大作图功能的软件系统,是由Ross Ihaka 和Robert Gentleman1共 同 创 立 。R语 言 可 以 看 作 是 由AT&T贝 尔 实 验 室 所 创 的S语 言 发 展 出 的 一 种 方 言1 。 因 此 ,R即 是 一 种 软 件 也 可 以 说 是 一 种 语 言 。S语 言 现 在 主 要 内 含 在 由Insightful2 公 司 经 营 的S-PLUS软 件 中 。R和S在 设计理念上存在有着许多不同:关于这方面的详细内容大家可以参考Ihaka & Gentleman (1996) 或R-FAQ3 ,该文档同时随R一起发布。 R是在GNU协议General Public Licence4 下免费发行的,它的开发及维护 现在则由R开发核心小组R Development Core Team 具体负责。 R的 安 装 文 件 有 多 种 形 式 , 有 在Unix 或Linux系 统 下 所 需 的 一 些 源 代 码(主要用C及Fortran 编写),及在Windows, Linux及Macintosh上使用的预编 译二进制码。这些安装文件以及安装说明都可以在Comprehensive R Archive Network (CRAN)5 网 站 上 下 载 。 该 网 站 提 供 的 关 于Linux的 安 装 文 件 只 适 用 于较新版本的Linux。详情请参考CRAN网站。 R内 含 了 许 多 实 用 的 统 计 分 析 及 作 图 函 数 。 作 图 函 数 能 将 产 生 的 图 片 展 示在一个独立的窗口中,并能将之保存为各种形式的文件(jpg, png, bmp, ps, pdf, emf, pictex, xfig; 具体形式取决于操作系统)。统计分析的结果也能被直 接 显 示 出 来 , 一 些 中 间 结 果(如P-值 , 回 归 系 数 , 残 差 等)既 可 保 存 到 专 门 的 文件中,也可以直接用作进一步的分析。 在R语言中,使用者可以使用循环语句来连续分析多个数据集,也可将多 个不同的统计函数结合在一个语句中执行更复杂的分析。R使用者还可以借鉴 网上提供的用S编写的大量程序6 ,而且大多数都能被R直接调用。 非专业人员起初可能觉得R相对比较复杂。其实,R的一个非常突出的优 点正是它的灵活性。一般的软件往往会直接展示分析的结果,而R则将这些结 果都存在一个对象“object”里面,所以常常在分析执行结束后并不显示任何结 1 Ihaka R. & Gentleman R. 1996. R: a language for data analysis and graphics. Journal of Computational and Graphical Statistics 5: 299–314. 2 http://www.insightful.com/products/splus/default.asp 3 http://cran.r-project.org/doc/FAQ/R-FAQ.html 4 For more information: http://www.gnu.org/ 5 http://cran.r-project.org/ 6 For example: http://stat.cmu.edu/S/ 1 果。使用者可能会对此感到困惑,其实这样的特点是非常有用的,因为我们 可以选择的从结果中只抽出我们感兴趣的部分。例如,我们要运行20个回归 分析而只想比较其回归系数,在R中就可以选择只显示所有分析得出的回归系 数,这样结果仅仅占了一排,而用有些软件可能会一下打开20个窗口。而在 下面的章节中,我们会看到更多能展示R相比传统软件更为灵活优越的例子。 2 2 基本原理与概念 如 果R已 经 被 安 装 在 你 的 计 算 机 中 , 它 就 能 立 即 运 行 一 些 可 执 行 的 命 令 了 。R默 认 的 命 令 提 示 符 是‘>’, 它 表 示 正 在 等 待 输 入 命 令 。 如 在Windows系 统中打开Rgui.exe,就能直接运行下拉菜单中的一些操作命令(如在线帮助, 打文件. . . )。到这里,有些人可能会急着想知道更多的语句命令。其实,在 学 习 这 些 内 容 前 , 了 解 掌 握 一 些R的 基 本 工 作 原 理 是 非 常 有 必 要 的 。 这 正 是 本 章所要讲的主要内容。 本 章 首 先 简 要 描 述R的 工 作 原 理 。 在 第 二 节 中 , 我 将 介 绍 一 些 基 本 的 赋 值分配(“assign”) 的操作,如怎样产生对象(object),如何操作管理这些对象 等.最后简要介绍R中非常有用的在线帮助。 2.1 基本原理 因为R是一种编程语言,一些对编程不太熟悉的人可能会望而却步。这种 障碍其实是完全没有必要,首先,R是一种解释型语言,而不是编译语言,也 就意味着输入的命令能够直接被执行,而不需要像一些语言要首先构成一个 完整的程序形式(如C,Fortan, Pascal, . . . )。 第二,R的语法非常之简单和直观。例如,线性回归的命令lm(y ~ x) 表 示“以x为 自 变 量 ,y 为 反 应 量来 拟 合 一 个 线 性 模 型”。 合 法 的R函 数 总 是 带 有 圆括号的形式,即使括号内没有内容(如,ls())。如果直接输入函数名而不输 入 圆 括 号 ,R 则 会 自 动 显 示 该 函 数 的 一 些 具 体 内 容 。 在 本 手 册 中 除 在 部 分 文 字 已作出清楚的说明外,所有的函数后都接有圆括号以区别于对象(object)。 当R运 行 时 , 所 有 变 量 , 数 据 , 函 数 及 结 果 都 以 对 象(objects )的 形 式 存 在计算机的活动内存中,并冠有相应的名字代号。我们可以通过用一些运算 符(如算术,逻辑,比较等)和一些函数(其本身也是对象)来对这些对象进行操 作。运算操作非常简单,其细节将留在下章讨论(p. 26). 关于R中的函数可用 下面的图例来形象的描述: arguments −→ options −→ function ↑ default arguments =⇒result 上 图 中 的 参 量(argument)可 能 是 一 些 对 象(如 数 据 , 方 程 , 算 式. . . )。 有 些参量在函数里被预设为缺省值,用户则可按需对其作个别的修改。所以运 行一个R函数可能不需要设定任何参量,原因是所有的参量都可以被默认为缺 3 省值,当然也有可能该函数本身就不含任何参量。由于这里主要是讲述R的工 作原理,对R函数的介绍将不再展开,在后面的章节中我们会看到关于构建及 使用各种函数的详细内容(p. 68)。 在R中进行的所有操作都是针对存储在活动内存中的对象的, 因此就不涉 及 到 任 何 临 时 文 件夹 的 使 用(Fig. 1)。 对 数 据 , 结 果 或 图 表 的 输 入 与 输 出 都 是通过在对计算机硬盘中的文件读写而实现。用户通过输入一些命令调用函 数,分析得出的结果可以被直接显示在屏幕上,也可以被存入某个对象或被 写入硬盘(如图片对象)。因为产生的结果本身就是一种对象,所以它们也能被 视为数据并能像一般数据那样被处理分析。数据文件即可从本地磁盘读取也 可通过网络传输从远程服务器端获得。 keyboard mouse commands - functions and operators .../library/base/ /stast/ /graphics/ ...  library of functions ?  “data” objects   )  screen data files -  XXX 6 XXX z X ? “results” objects PS Active memory internet JPEG ... Hard disk 图 1: R工作原理示意图. 所 有 能 使 用 的R函 数 都 被 包 含 在 一 个 库(library) 中, 该 库 存 放 在 磁 盘 的R HOME/library 目 录 下(R HOME 是 最 初 安 装R的地 址)。 这 个 目 录 下 含 有 具 有 各 种 功 能 的 包(packages ), 各个 包 也 是 按 照 目 录 的 方 式 组 织 起 来 的 。 其中名为base的包可以算是R的核心,因为它内嵌了R语言中所有像数据读写 与 操 作 这 些 最 基 本 的 函 数 。 在 上 述 目 录 中 的 每 个 包 内 , 都 有 一 个 子 目 录R, 这个目录里又都含有一个与包同名的文件(例如在包base中,有这样一个文 件R HOME/library/base/R/base)。该文件正是存放所有函数的地方。 R语 言 中 最 简 单 的 命 令 莫 过 于 通 过 输 入 一 个 对 象 的 名 字 来 显 示 其 内 容 了 。 例如,一个名为n的对象,其内容是数值10: >n [1] 10 方括号中的数字1表示从n的第一个元素开始显示。其实该命令的功能在 4 这里于函数print 相似,输出结果与print(n) 相同(但有些情况下,例如内嵌 在一个函数或循环中时,就必须得用print函数)。 对象的名字必须是以一个字母开头(A–Z 或a–z), 中间可以包含字母,数 字(0–9),点(.)及下划线( ). 因为R对对象的名字区分大小写,所以x 和X就可 以代表两个完全不同的对象(在Windows操作系统中也是如此)。 对象的产生,排列及删除 2.2 一个对象可以通过赋值操作来产生,R语言中的赋值(“assign”) 符号一般 是由一个尖括号与一个负号组成的箭头形标志。该符号可以是从左到右的方 向,也可以相反: >n >n [1] >5 >n [1] >x >X >x [1] >X [1] <- 15 15 -> n 5 <- 1 <- 10 1 10 如果该对象已经存在,那么它以前的值将会自动被新值冲掉(这种修改只 会影响内存中的数据,操作结果暂时不会被保存到硬盘中)。在R中给对象赋 值有多种形式,可以是直接赋一个数值,也可以是一个算式或一个函数的结 果: >n >n [1] >n >n [1] <- 10 + 2 12 <- 3 + rnorm(1) 2.208807 运行rnorm(1) 将产生一个服从平均数为0标准差为1的标准正态分布的随 机变量(p. 18)。当然你也可以只是输入函数或表达式而不把它的结果赋给某 个对象,但如果这样在窗口中展示的结果将不会被保存到内存中: 5 > (10 + 2) * 5 [1] 60 本文中,在不影响读者理解的情况下,一些赋值符号将会被省略。 函数ls的功能是显示所有在内存中的对象:只会列出对象名,如: > name <- "Carmen"; n1 <- 10; n2 <- 100; m <- 0.5 > ls() [1] "m" "n1" "n2" "name" 注 意 在R 中 应 该 用 分 号 来 隔 开 同 一 行 中 的 不 同 命 令 语 句 。 如 果 只 要 显 示 出 在名称中带有某个指定字符的对象,则通过设定选项pattern 来实现(可简写 为pat) ): > ls(pat = "m") [1] "m" "name" 如果进一步限为显示在名称中以某个字母开头的对象,则可: > ls(pat = "^m") [1] "m" 运行函数ls.str()将会展示内存中所有对象的详细信息: > ls.str() m : num 0.5 n1 : num 10 n2 : num 100 name : chr "Carmen" 选 项pattern在这 里 同 样 适 用 。 在ls.str函 数 中 另 一 个 非 常 有用 的 选 项 是max.level, 它 将 规 定 显 示 所 有 关 对 象 信 息 的 详 细 级 别 。 缺 省 情 况 下 ,ls.str 将 会 列 出 关 于 对 象 的 所 有 信 息 , 包 括 数 据 框 , 矩 阵 , 数 据 列 表的列数信息。因此展示结果可能会很长。但如果设定max.level =-1 就可 以避免这种情况了: > M <- data.frame(n1, n2, m) > ls.str(pat = "M") M : ‘data.frame’: 1 obs. of $ n1: num 10 $ n2: num 100 $ m : num 0.5 > ls.str(pat="M", max.level=-1) M : ‘data.frame’: 1 obs. of 3 variables: 3 variables: 要在内存中删除某个对象,可利用函数rm: 运行rm(x)将会删除对象x,运 行rm(x,y) 将会删除对象x和y,而运行rm(list=ls())则会删除内存中的所有对 象. 在上面所讲的ls() 函数中的一些选项同样也可以运用到rm中来,以选择 的删除一些对象,如: rm(list=ls(pat="^m"))。 6 2.3 在线帮助 R中给予的在线帮助能提供关于如何使用函数的非常有用的信息。关于某 个特定函数的帮助能够直接被调出来,如运行: > ?lm 会立即显示关于函数lm()(线性模型)的帮助页面。命令help(lm) 和help("lm") 具有同样的效果。但在查询关于某特殊语法意义字符的帮助时必须用后一种 形式,如: > ?* Error: syntax error > help("*") Arithmetic package:base R Documentation Arithmetic Operators ... 启动帮助将会打开一个页面(取决于操作系统),第一行一般会显示某函数 或操作命令的所属的包(package),然后是标题,标题下面是则是一些详细信 息。 Description: brief description. Usage: for a function, gives the name with all its arguments and the possible options (with the corresponding default values); for an operator gives the typical use. Arguments: for a function, details each of its arguments. Details: detailed description. Value: if applicable, the type of object returned by the function or the operator. See Also: other help pages close or similar to the present one. Examples: some examples which can generally be executed without opening the help with the function example. 对 初 学 者 而 言 , 参 考 帮 助中Examples部 分 的 信 息 是 很 有用 的 。 而 一 般 应该仔细阅读Arguments中的一些说明也是非常有必要的。帮助中还包含了 其它一些说明部分,如Note, References或Author(s)等。 默认状态下,函数help只会在被载入内存中的包中搜索。选项try.all.package 在缺省值是FALSE, 但如果把它设为TRUE,则可在所有包中进行搜索: 7 > help("bs") No documentation for ’bs’ in specified packages and libraries: you could try ’help.search("bs")’ > help("bs", try.all.packages = TRUE) Help for topic ’bs’ is not in any loaded package but can be found in the following packages: Package splines Library /usr/lib/R/library 但注意在这种情况下,不会显示关于函数bs的帮助页面,如果使用者确 实想打开这样的页面而所属包又没有被载入内存时,可以使用package这个选 项: > help("bs", package = "splines") bs package:splines R Documentation B-Spline Basis for Polynomial Splines Description: Generate the B-spline basis matrix for a polynomial spline. ... Html格式的帮助可以通过输入下面的函数启动 > help.start() 在html格式的帮助页面中还可以使用关键词进行搜索。在See Also部分 中 , 可 以 通 过 超 文 本 链 接 到 其 他 相 关函 数 的 帮 助 页 面 。 使 用 关 键 词 的 搜索 在R中 也 可 以 通 过函 数help.search来 实 现 。 这 种 方法 能 在 所 有 已 安 装 的 包 中 搜索 包 含 给 定 字 符 串 的 相 关 内 容 。 例 如 , 运 行help.search("tree")会 列 出 所 有 在 帮 助 页 面 含 有“tree”的 函 数 。 注 意 如 果 有 一 些 包 是 最 近 才 安 装 的 , 应 该 首 先 使 用 函 数help.search中 的rebuild选 项 来 刷 新 数 据 库(e.g., help.search("tree", rebuild = TRUE))。 使用函数apropos则能找出所有在名字中含有指定字符串的函数,但只会 在被载入内存中的包中进行搜索: > apropos(help) [1] "help" [4] "help.start" ".helpForCall" "help.search" 8 3 3.1 R的数据 操作 对象 我们已经看到R通过一些对象来运行,当然首先这些对象是用它们的名称 和内容来刻画的,其次也通过对象的数据类型即属性来刻画。为了理解这些 属性的用处,我们以一个在{1,2,3}中取值的变量为例:这个变量可以是一个 整数变量(例如巢中蛋的个数),或者也可以是一个分类变量的编码(例如 某些甲壳类动物的三种性别:雄、雌和雌雄同体)。 显 然 对 这 个 变 量 的 统 计 分 析 在 以 上 两例 中 将 是 不 相 同 的 , 对 象 的 属 性 在R中提供着所需的信息。更技术性也更一般地说,对于作用于一个对象的函 数,其表现将取决于对象的属性。 所有的对象都有两个内在属性:类型和长度。类型是对象元素的基本种 类 , 共 有 四 种 : 数 值 型 , 字 符 型 , 复 数 型7 和 逻 辑 型(FALSE或TRUE), 虽 然 也 存在其它的类型,但是并不能用来表示数据,例如函数或表达式;长度是对 象中元素的数目。对象的类型和长度可以分别通过函数mode和length得到。 > x <- 1 > mode(x) [1] "numeric" > length(x) [1] 1 > A <- "Gomphotherium"; compar <- TRUE; z <- 1i > mode(A); mode(compar); mode(z) [1] "character" [1] "logical" [1] "complex" 无 论 什 么 类 型 的 数 据 , 缺 失数 据 总 是 用NA(不可用)来 表 示 ; 对 很 大 的 数 值则可用指数形式表示: > N <- 2.1e23 >N [1] 2.1e+23 R可 以 正 确 地 表 示 无 穷 的 数 值 , 如 用Inf和-Inf表 示±∞, 或 者 用NaN(非 数字)表示不是数字的值。 7 本手册中不讨论复数型 9 > x <- 5/0 >x [1] Inf > exp(x) [1] Inf > exp(-x) [1] 0 >x-x [1] NaN 字 符 型 的 值 输 入 时 须 加 上 双 引 号", 如 果 需 要引 用 双 引 号 的 话 , 可 以 让 它 跟 在 反 斜 杠\后 面 ; 这 两 个 字 符 合 一 起\"在 某 些 函 数 如cat的 输 出 显 示 或write.table写入磁盘(参见p. 14,函数的qmethod选项)时会被以特殊的方式 处理。 > x <- "Double quotes \" delimitate R’s strings." >x [1] "Double quotes \" delimitate R’s strings." > cat(x) Double quotes " delimitate R’s strings. 也有另一种表示字符型变量的方法,即用单引号(’)来界定变量,这种情 况下不需要用反斜杠来引用双引号(但是引用单引号时必须要用!) > x <- ’Double quotes " delimitate R\’s strings.’ >x [1] "Double quotes \" delimitate R’s strings." 下表给出了表示数据的对象的类别概览: 对象 向量 因子 数组 矩阵 数据框 时间序列(ts) 列表 类型 是否允许 同一个对象中 有多种类型? 数值型,字符型,复数型,或 逻辑型 数值型或 字符型 数值型,字符型,复数型,或 逻辑型 数值型,字符型,复数型,或 逻辑型 数值型,字符型,复数型,或 逻辑型 数值型,字符型,复数型,或 逻辑型 数值型,字符型,复数型,逻辑型, 函数,表达式,. . . 否 否 否 否 是 否 是 10 向量是一个变量,其意思也即人们通常认为的那样;因子是一个分类变 量 ; 数 组 是 一 个k 维 的 数 据 表 ; 矩 阵 是数 组 的 一 个 特 例 , 其 维 数k = 2。 注 意 , 数 组 或 者 矩 阵中 的 所 有 元 素 都 必 须 是 同 一 种 类 型 的 ; 数 据 框 是 由 一 个 或几个向量和(或)因子构成,它们必须是等长的,但可以是不同的数据类 型 ; “ts” 表 示时 间 序 列 数 据 , 它 包 含 一 些 额 外 的 属 性 , 例 如 频 率 和 时 间 ; 列表可以包含任何类型的对象,包括列表! 对于一个向量,用它的类型和长度足够描述数据;而对其它的对象则另 需一些额外信息,这些信息由外在的属性给出。这些属性中的是表示对象维 数 的dim , 比 如 一 个2行2列 的的 矩 阵 , 它 的dim是 一 对 数 值[2,2], 但 是 其 长 度 是4。 3.2 在文件中读写数据 对 于 在 文 件 读 取 和 写 入 的 工 作 ,R使 用 工 作 目 录 来 完 成 。 可 以 使 用 命 令getwd() (获 得 工 作 目 录 )来 找 到 目 录 , 使 用 命 令setwd("C:/data") 或 者setwd("/home/paradis/R") 来 改 变 目 录 。 如 果 一 个 文 件 不 在 工 作 目 录 里 则必须给出它的路径8 。 R可以用下面的函数读取存储在文本文件(ASCII)中的数据:read.table (其 中 有 若 干 参 数 , 见 后 文),scan和read.fwf。R也 可 以 读 取 以 其 他 格 式 的 文 件(Excel, SAS, SPSS, . . . ) 和 访 问SQL类 型 的 数 据 库 , 但 是 基 础 包 中 并不 包含所需的这些函数。这些功能函数对于R的高级应用是十分有用的,但是我 们在这里将读取文件限定在ASCII格式。 函数read.table用来创建一个数据框,所以它是读取表格形式的数据的 主要方法。举例来说,对于一个名为data.dat的文件,命令: > mydata <- read.table("data.dat") 将 创 建 一 个 数 据 框 名 为mydata, 数 据 框 中 每 个 变 量 也 都 将 被 命 名 , 缺 省 值 为V1, V2, . . . 并 且 可 以 单 独地 访 问 每 个 变 量 代 码 为 :mydata$V1, mydata$V2, . . . , 或 者 用mydata["V1"], mydata["V2"], . . . , 或 者 还 有 一 种 方 法,mydata[, 1], mydata[,2 ], . . . 9 这里有一些选项的缺省值(即如果用户 不设定那么R将自动使用的值)见于下表: read.table(file, header = FALSE, sep = "", quote = "\"’", dec = ".", row.names, col.names, as.is = FALSE, na.strings = "NA", colClasses = NA, nrows = -1, skip = 0, check.names = TRUE, fill = !blank.lines.skip, 8 在Windows中,为Rgui.exe创建一个快捷方式是比较有用的,在快捷方式“属性”的“起 始位置”中改变目录,然后用此快捷方式启动R时这个目录就会成为工作目录 9 注 意 这 几 种 方 法 的 结 果 是 有 区 别 的 :mydata$V1和mydata[, 1]是 向 量 , 而mydata["V1"]是数据框。后面(p. 19)将会讲到关于处理对象的详情。 11 strip.white = FALSE, blank.lines.skip = TRUE, comment.char = "#") file header sep quote dec row.names col.names as.is na.strings colClasses nrows skip check.names fill strip.white blank.lines.skip comment.char 文 件 名 ( 包 在""内 , 或 使 用 一 个 字 符 型 变 量 ), 可 能 需 要 全 路 径 (注意即使是在Windows下,符号\ 也不允许包含在内,必须用/替 换),或者一个URL链接(http://...)(用URL对文件远程访问) 一 个 逻 辑 值(FALSE or TRUE), 用 来 反 映 这 个 文 件 的第 一 行 是 否 包 含 变量名 文件中的字段分离符,例如对用制表符分隔的文件使用sep="\t" 指定用于包围字符型数据的字符 用来表示小数点的字符 保 存 着 行 名 的 向 量,或 文 件 中 一 个 变 量 的 序 号 或 名 字,缺 省时 行 号 取 为1, 2, 3, . . . 指定列名的字符型向量(缺省值是:V1, V2, V3, . . . ) 控制是否将字符型变量转化为因子型变量(如果值为FALSE),或者仍 将其保留为字符型(TRUE)。as.is可以是逻辑型,数值型或者字符 型向量,用来判断变量是否被保留为字符。 代表缺失数据的值(转化为NA) 指定各列的数据类型的一个字符型向量 可以读取的最大行数(忽略负值) 在读取数据前跳过的行数 如果为TRUE,则检查变量名是否在R中有效 如果为TRUE且非所有的行中变量数目相同,则用空白填补 在sep已指定的情况下,如果为TRUE,则删除字符型变量前后多余的 空格 如果为TRUE,忽略空白行 一个字符用来在数据文件中写注释,以这个字符开头的行将被忽略 (要禁用这个参数,可使用comment.char = "") read.table的几个变种因为使用了不同的缺省值可以用在几种不同情况 下: read.csv(file, header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, ...) read.csv2(file, header = TRUE, sep = ";", quote="\"", dec=",", fill = TRUE, ...) read.delim(file, header = TRUE, sep = "\t", quote="\"", dec=".", fill = TRUE, ...) read.delim2(file, header = TRUE, sep = "\t", quote="\"", dec=",", fill = TRUE, ...) 函数scan比read.table要更加灵活,它们的区别之一是前者可以指定变 量的类型,例如: 12 > mydata <- scan("data.dat", what = list("", 0, 0)) 读 取 了 文 件data.dat中 三 个 变 量 , 第 一 个 是 字 符 型 变 量 , 后 两 个 是数 值 型变量。另一个重要的区别在于scan()可以用来创建不同的对象,向量,矩 阵,数据框,列表. . . 在上面的例子中,mydata是一个有三个向量的列表。 在 缺 省 情 况 下 , 也 就 是 说 , 如 果what被 省 略 ,scan()将 创 建 一 个 数 值 型 向 量。如果读取的数据类型与缺省类型或指定类型不符,则将返回一个错误信 息。这些选项在下面进行说明。 scan(file = "", what = double(0), nmax = -1, n = -1, sep = "", quote = if (sep=="\n") "" else "’\"", dec = ".", skip = 0, nlines = 0, na.strings = "NA", flush = FALSE, fill = FALSE, strip.white = FALSE, quiet = FALSE, blank.lines.skip = TRUE, multi.line = TRUE, comment.char = "") file what nmax n sep quote dec skip nlines na.string flush fill strip.white quiet blank.lines.skip multi.line comment.char 文 件 名(在""之 内), 可 能 包 含 它 的 路 径 ( 符 号\ 不 允 许 使 用 , 必 须 用/替 代, 即 使 是 在Windows下 面 ) , 或 者 使 用 一 个URL链 接 (http://...) ; 如 果file=””, 数 据 从 键 盘 输 入 ( 使 用 一 个 空 白 行 终止输入) 指定数据的类型(缺省值为数值型) 要 读 取 数 据 的 最 大 数 量 , 如 果what是 一 个 列 表 ,nmax则 是 可 以 读 取 的行数(在缺省情况下,scan读取到文件最末端为止的所有数据) 要读取数据的最大数量(在缺省情况下,没有限制) 文件中的字段分隔符 用来包围字符型值 用来表示小数点的字符 在读取数据前跳过的行数 要读取的行数 表示缺失数据的字符串(转化为为NA) 一 个 逻 辑 值 , 如 果 为TRUE, 当读 取 完 指 定 列 数 后scan将 转 到 下 一 行 (这样就允许用户在数据文件中添加注释,即添加在指定列数之后) 如果为TRUE,且非所有的行中变量数目相同,则用空白填补 在sep已指定的情况下,如果为TRUE,则删除字符型变量前后多余的 空格 一个逻辑值,如果为FALSE,scan显示一行信息说明哪些字段被读取 如果为TRUE,忽略空白行 当what是一个列表时,若为FALSE则表示列表中每个个体的所有变量 都在同一行中 指定注释开始字符,一行中以这个字符开头的部分将被忽略(缺是禁 用此选项) 函数read.fwf可以用来读取文件中一些固定宽度格式的数据: read.fwf(file, widths, sep="\t", as.is = FALSE, skip = 0, row.names, col.names, n = -1, ...) 13 除 了widths用 来 说 明 读 取 字 段 的 宽 度 外 , 选 项 与read.table()基 本 相 同 。 举 例来 说 , 如 果 在 一 个 名 为data.txt的文件中有一组如右所示的数据,可以读取这 些数据用下面的命令: A1.501.2 A1.551.3 B1.601.4 B1.651.5 C1.701.6 C1.751.7 > mydata <- read.fwf("data.txt", widths=c(1, 4, 3)) > mydata V1 V2 V3 1 A 1.50 1.2 2 A 1.55 1.3 3 B 1.60 1.4 4 B 1.65 1.5 5 C 1.70 1.6 6 C 1.75 1.7 3.3 存储数据 函数write.table可以在文件中写入一个对象,一般是写一个数据框,也 可以是其它类型的对象(向量,矩阵. . . )。参数和选项: write.table(x, file = "", append = FALSE, quote = TRUE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = TRUE, col.names = TRUE, qmethod = c("escape", "double")) 14 x file append quote sep eol na dec row.names col.names qmethod 要写入的对象的名称 文件名(缺省时对象直接被“写”在屏幕上) 如果为TRUE则在写入数据时不删除目标文件中可能已存在的数据,采取往后 添加的方式 一 个 逻 辑 型 或 者 数 值 型 向 量 : 如 果 为TRUE, 则 字 符 型 变 量 和 因 子 写 在 双 引 号""中;若quote是数值型向量则代表将欲写在""中的那些列的列标。(两种 情况下变量名都会被写在""中;若quote = FALSE则变量名不包含在双引号 中) 文件中的字段分隔符 使用在每行最后的字符("\n"表示回车) 表示缺失数据的字符 用来表示小数点的字符 一个逻辑值,决定行名是否写入文件;或指定要作为行名写入文件的字符型 向量 一个逻辑值(决定列名是否写入文件);或指定一个要作为列名写入文件中 的字符型向量 若quote=TRUE, 则 此 参 数 用 来 指 定 字 符 型 变 量 中 的 双 引 号"如 何 处 理 : 若 参 数 值 为"escape" (或 者"e", 缺 省)每 个"都 用\"替 换 ; 若 值 为"d"则 每 个"用""替换 若 想 用 更 简 单 的 方 法 将 一 个 对 象 写 入 文 件 , 可 以 使 用 命 令write(x, file ="data.txt"),其 中x是 对 象 的 名 字 ( 它 可 以 是 向 量 , 矩 阵 , 或 者 数 组 ) 。 这 里 有 两 个 选 项 :nc(或 者ncol), 用 来 定 义 文 件 中 的 列 数 ( 在 缺 省 情 况 下 , 如 果x是 字 符 型 数 据 , 则nc=1; 对 于 其 它 数 据 类 型nc=5) , 和append( 一 个 逻 辑 值 ), 若 为TRUE则 添 加 数 据 时 不 删 除 那 些 可 能 已 存 在 在文件中的数据;若为FALSE(缺省值)则删除文件中已存在的数据。 要 记 录 一 组 任 意 数 据 类 型 的对 象 , 我 们 可 以 使 用 命 令save(x, y, z, file= "xyz.RData")。 可 以 使 用 选 项ASCII=TRUE使 得 数 据 在 不 同 的 机 器 之 间更简易转移。数据(用R的术语来说叫做工作空间)可以在使用load("xyz. RData")之后被加载到内存中。函数save.image()是save(list =ls(all=TRUE), file=".RData")的一个简捷方式。 生成数据 3.4 3.4.1 规则序列 例如一个从1到30的规则整数序列,可以这样产生: > x <- 1:30 这个结果向量x有30个元素。算子‘:’的优先级可从如下表达式中看出: > 1:10-1 [1] 0 1 2 3 4 5 6 7 8 9 > 1:(10-1) 15 [1] 1 2 3 4 5 6 7 8 9 函数seq可以生成如下的实数序列: > seq(1, 5, 0.5) [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 其中第一个数字表示序列的起点,第二个表示终点,第三个是生成序列的步 长。也可以这样使用: > seq(length=9, from=1, to=5) [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 还可以用函数c直接输入数值: > c(1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5) [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 如果想用键盘输入一些数据也是可以的,只需要直接使用默认选项 的scan函数: > z <- scan() 1: 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 10: Read 9 items >z [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 函数rep用来创建一个所有元素都相同的向量: > rep(1, 30) [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 函数sequence创建一系列连续的整数序列,每个序列都以给定参数的数 值结尾: > sequence(4:5) [1] 1 2 3 4 1 2 3 4 5 > sequence(c(10,5)) [1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 函 数gl(生成不同的水平/层次数据)十 分 有用 , 因 为 它 能 产 生 规 则 的 因 子 序 列 。 这 个 函 数 的 用 法 是gl(k,n), 其 中k是 水 平 数 ( 或 类 别 数 ),n是 每 个 水 平 重 复 的 次 数 。 此 函 数 有 两 个 选 项 :length用 来 指 定 产 生数 据 的 个 数,labels用来指定每个水平因子的名字。例如: 16 > gl(3, 5) [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 Levels: 1 2 3 > gl(3, 5, length=30) [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 Levels: 1 2 3 > gl(2, 6, label=c("Male", "Female")) [1] Male Male Male Male Male Male [7] Female Female Female Female Female Female Levels: Male Female > gl(2, 10) [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 Levels: 1 2 > gl(2, 1, length=20) [1] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 Levels: 1 2 > gl(2, 2, length=20) [1] 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 Levels: 1 2 最 后 ,expand.grid()创 建 一 个 数 据 框 , 结 果 是 把 各 参 数 的 各 水 平 完 全 搭配: > expand.grid(h=c(60,80), w=c(100, 300), sex=c("Male", "Female")) h w sex 1 60 100 Male 2 80 100 Male 3 60 300 Male 4 80 300 Male 5 60 100 Female 6 80 100 Female 7 60 300 Female 8 80 300 Female 17 3.4.2 随机序列 分布名称 函数 Gaussian (normal) exponential gamma Poisson Weibull Cauchy beta ‘Student’ (t) Fisher–Snedecor (F ) Pearson (χ2 ) binomial multinomial geometric hypergeometric logistic lognormal negative binomial uniform Wilcoxon’s statistics rnorm(n, mean=0, sd=1) rexp(n, rate=1) rgamma(n, shape, scale=1) rpois(n, lambda) rweibull(n, shape, scale=1) rcauchy(n, location=0, scale=1) rbeta(n, shape1, shape2) rt(n, df) rf(n, df1, df2) rchisq(n, df) rbinom(n, size, prob) rmultinom(n, size, prob) rgeom(n, prob) rhyper(nn, m, n, k) rlogis(n, location=0, scale=1) rlnorm(n, meanlog=0, sdlog=1) rnbinom(n, size, prob) runif(n, min=0, max=1) rwilcox(nn, m, n), rsignrank(nn, n) 在统计学中,产生随机数据是很有用的,R可以产生多种不同分布下的随 机 数 序 列 。 这 些 分 布 函 数 的 形 式 为rfunc(n,p1,p2,...), 其 中func指 概 率 分 布函数,n为生成数据的个数,p1, p2, . . . 是分布的参数数值。上面的表给出 了每个分布的详情和可能的缺省值(如果没有给出缺省值,则意味着用户必 须指定参数)。 大多数这种统计函数都有相似的形式,只需用d、p或者q去替代r,比如 密度函数(dfunc (x, ...)),累计概率密度函数(也即分布函数)(pfunc (x, ...))和分位数函数(qfunc (p, ...),0 < p < 1)。最后两个函数序列可以用 来求统计假设检验中P值或临界值。例如,显著性水平为5%的正态分布的双 侧临界值是: > qnorm(0.025) [1] -1.959964 > qnorm(0.975) [1] 1.959964 对于同一个检验的单侧临界值,根据备择假设的形式使用qnorm(0.05)或1 qnorm(0.95)。 18 一个检验的P 值,比如自由度df = 1的χ2 = 3.84: > 1 - pchisq(3.84, 1) [1] 0.05004352 3.5 3.5.1 使用对象 创建对象 我们在前面看到了用赋值操作创建对象的不同方法;在这样的创建中对 象的数据类型和模式通常都已经预先确定了。在创建一个对象时是有可能指 定它的数据类型、长度、类别等等的。从处理对象的角度来看这些方法是很 有趣的。举例来说,我们可以创建一个空的对象并且逐步修改其中的元素, 这比把所有的元素一起用c()放进去更有效。在这里也可以使用下标系统,后 面我们将会看到(p. 27)。 从一些对象创建其它对象也是很方便的,比如,若想拟合一系列线性模 型,可以将这一系列模型的公式放在一个列表中,然后不断地取出元素(即 那些公式)插入到函数lm中。 在这个R的学习阶段,我们很有必要了解下面的函数和数据结构。直接创 建数据结构不仅能让我们对数据有更好的理解,而且也会更深入地领会前文 中提到的一些概念。 向量( Vector) 函 数vector有 两 个 参 数 : 类 型(mode)和 长 度(length), 创 建的向量中元素值取决于参数所指定的数据类型:数值型向量则元素值 都为0,逻辑型都为FALSE,字符型都为""。以下三个函数有几乎相同的 效果(创建一个向量)并且只有一个参数即长度:numeric(),logical(), 和character()。 因子( Factor) 一 个 因 子 不 仅 包 括 分 类 变 量 本 身 还 包 括 变 量 不 同 的 可 能 水 平(即使它们在数据中不出现)。因子函数factor用下面的选项创建一 个因子: factor(x, levels = sort(unique(x), na.last = TRUE), labels = levels, exclude = NA, ordered = is.ordered(x)) levels 用来指定因子可能的水平(缺省值是向量x中互异的值);labels 用来指定水平的名字;exclude表示从向量x中剔除的水平值;ordered是 一个逻辑型选项用来指定因子的水平是否有次序。回想数值型或字符型 的x。下面有一些例子: > factor(1:3) [1] 1 2 3 19 Levels: 1 2 3 > factor(1:3, levels=1:5) [1] 1 2 3 Levels: 1 2 3 4 5 > factor(1:3, labels=c("A", "B", "C")) [1] A B C Levels: A B C > factor(1:5, exclude=4) [1] 1 2 3 NA 5 Levels: 1 2 3 5 函数levels用来提取一个因子中可能的水平值: > ff <- factor(c(2, 4), levels=2:5) > ff [1] 2 4 Levels: 2 3 4 5 > levels(ff) [1] "2" "3" "4" "5" 矩阵( Matrix) 一 个 矩 阵 实 际 上 是 有 一 个附 加 属 性 ( 维 数dim) 的 向 量 , 维数即为一个长度为2的向量,用来指定矩阵的行数和列数。一个矩阵 可以用函数matrix来创建: matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL) 选 项byrow表 示数 据 给 出 的 值 是 要 按 列 填 充 ( 缺 省 值 ) 还 是 按 行 填 充 (如果为TRUE)。可以通过选项dimnames给行列命名。 > matrix(data=5, nr=2, nc=2) [,1] [,2] [1,] 5 5 [2,] 5 5 > matrix(1:6, 2, 3) [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 > matrix(1:6, 2, 3, byrow=TRUE) [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 20 另一种创建矩阵的方法是给维数适当赋值(初始值为NULL): > x <- 1:15 >x [1] 1 2 3 4 5 > dim(x) NULL > dim(x) <- c(5, 3) >x [,1] [,2] [,3] [1,] 1 6 11 [2,] 2 7 12 [3,] 3 8 13 [4,] 4 9 14 [5,] 5 10 15 6 7 8 9 10 11 12 13 14 15 数据框 (Data frame) 前 面 我 们 已 经 看 到 一 个 数 据 框 可 以 由 函 数read. table 间接创建;这里也可以用函数data.frame来创建。数据框中的向 量必须有相同的长度,如果其中有一个比其它的短,它将“循环”整数 次(以使得其长度与其它向量相同): > x <- 1:4; n <- 10; M <- c(10, 35); y <- 2:4 > data.frame(x, n) xn 1 1 10 2 2 10 3 3 10 4 4 10 > data.frame(x, M) xM 1 1 10 2 2 35 3 3 10 4 4 35 > data.frame(x, y) Error in data.frame(x, y) : arguments imply differing number of rows: 4, 3 如果一个因子包含在一个数据框中,它必须和其中的向量有相同的长 度 。 列 名 也 是 可 以 改 变 的 , 例 如 ,data.frame(A1=x, A2=n)。 用 户 也 21 可以使用row.names给行命名,但是,这个命名向量必须是字符型的而 且长度等于这个数据框的行数。最后,注意数据框和矩阵一样有维数这 个属性。 列表( List) 列表可以用list函数创建,方法与创建数据框类似。它对其中 包 含 的对 象 没 有 什 么 限 制 。 和data.frame()比 较 , 缺 省 值 没 有 给 出 对 象的名称;用前面的向量x和y举例: > L1 <- list(x, y); L2 <- list(A=x, B=y) > L1 [[1]] [1] 1 2 3 4 [[2]] [1] 2 3 4 > L2 $A [1] 1 2 3 4 $B [1] 2 3 4 > names(L1) NULL > names(L2) [1] "A" "B" 时间序 列(Time-series) 函 数ts可 以 由 向 量 ( 一 元 时 间 序 列 ) 或 者 矩 阵 (多元时间序列)创建一个ts型对象,并且有一些表明序列特征的选项 (带有缺省值),它们是: ts(data = NA, start = 1, end = numeric(0), frequency = 1, deltat = 1, ts.eps = getOption("ts.eps"), class, names) 22 data start end frequency deltat ts.eps class names 一个向量或者矩阵 第一个观察值的时间,为一个数字或者是一个由两个整 数构成的向量(参见下面的例子) 最后一个观察值的时间,指定方法和start相同 单位时间内观察值的频数(频率) 两个观察值间的时间间隔(例如,月度数据的取值 为1/12);frequency和deltat必 须 并 且 只 能 给 定 其 中 之一 序列之间的误差限。如果序列之间的频率差异小 于ts.eps则认为这些序列的频率相等。 对象的类型;一元序列的缺省值是"ts",多元序列的缺 省值是c("mts", "ts") 一个字符型向量,给出多元序列中每个一元序列的 名 称 ; 缺 省 为data中 每 列 数 据 的 名 称 或 者Series 1, Series 2, . . . 用ts创建时间序列的一些例子: > ts(1:10, start = 1959) Time Series: Start = 1959 End = 1968 Frequency = 1 [1] 1 2 3 4 5 6 7 8 9 10 > ts(1:47, frequency = 12, start = c(1959, 2)) Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 1959 1 2 3 4 5 6 7 8 9 10 11 1960 12 13 14 15 16 17 18 19 20 21 22 23 1961 24 25 26 27 28 29 30 31 32 33 34 35 1962 36 37 38 39 40 41 42 43 44 45 46 47 > ts(1:10, frequency = 4, start = c(1959, 2)) Qtr1 Qtr2 Qtr3 Qtr4 1959 1 2 3 1960 4 5 6 7 1961 8 9 10 > ts(matrix(rpois(36, 5), 12, 3), start=c(1961, 1), frequency=12) Series 1 Series 2 Series 3 Jan 1961 8 5 4 Feb 1961 6 6 9 Mar 1961 2 3 3 23 Apr May Jun Jul Aug Sep Oct Nov Dec 1961 1961 1961 1961 1961 1961 1961 1961 1961 8 4 4 4 11 6 6 5 8 5 9 6 2 6 5 5 5 5 4 3 13 6 4 7 7 7 2 表达式 (Expression) 表达式类型的对象在R中有着很基础的地位,是R能 够解释的字符序列。所有有效的命令都是表达式。一个命令被直接从键 盘输入后,它将被R求值,如果是有效的则会被执行。在很多情况下, 构造一个不被求值的表达式是很有用的:这就是函数expresssion要做 的。当然也可以随后用eval()对创建的表达式进行求值。 > x <- 3; y <- 2.5; z <- 1 > exp1 <- expression(x / (y + exp(z))) > exp1 expression(x/(y + exp(z))) > eval(exp1) [1] 0.5749019 达式也可以在其它地方用来在图表中添加公式(p. 43);表达式可以由 字符型变量创建;一些函数把表达式当作参数,例如可以求偏导数的函 数D。 > D(exp1, "x") 1/(y + exp(z)) > D(exp1, "y") -x/(y + exp(z))^2 > D(exp1, "z") -x * exp(z)/(y + exp(z))^2 3.5.2 对象的类型转换 读者必然会发现一些类型的对象之间的差异是很小的;因此改变一个对 象的某些属性使它转换为另一种类型的对象是合乎逻辑的。as.something 这 种形式的函数可以完成转换。R(2.1.0版)的base和utils包中有98个这种函数 在里面,所以我们在这里不做深入的阐述。 24 很明显转换取决于被转换对象的属性。一般来说,转换遵循一些很直观 的规则。对于类型的不同转换,下表总结了不同的情况。 转换目标 函数 数值型 as.numeric 逻辑型 as.logical 字符型 as.character 规则 FALSE TRUE "1", "2", . . . "A", . . . 0 其它数字 "FALSE", "F" "TRUE", "T" 其它字符 1, 2, . . . FALSE TRUE → → → → → → → → → → → → 0 1 1, 2, . . . NA FALSE TRUE FALSE TRUE NA "1", "2", . . . "FALSE" "TRUE" 有许多函数可以用来转换对象的类型(as.matrix, as.ts, as.data.frame, as.expression, . . . ),这些函数在转换时会影响除了类型之外的属性,将得 到的结果在一般情况下也是容易预见的。将因子转换为数值型是R中经常遇到 的情况,这种情况下R将因子的水平转化为数值编码: > fac <- factor(c(1, 10)) > fac [1] 1 10 Levels: 1 10 > as.numeric(fac) [1] 1 2 这对于一个字符型因子来说是很有意义的: > fac2 <- factor(c("Male", "Female")) > fac2 [1] Male Female Levels: Female Male > as.numeric(fac2) [1] 2 1 注意这个结果可能不像根据上面表格预期的那样是NA。 要想将一个数值型因子转换为一个数值型向量并且保持最初指定的水平 值,就必须先转换成字符型然后再转换成数值型。 25 > as.numeric(as.character(fac)) [1] 1 10 当一个文件中的数值型变量具有非数值的值时,这种做法就很有用;前 面我们已经看到了read.table()在缺省情况下会读取这列为一个因子。 3.5.3 运算符 我们在前面已经看到了R10 中三种主要类型的运算符。下面是运算符的列 表: 运算符 比较运算 数学运算 + * / ^ %% %/% 加法 减法 乘法 除法 乘方 模 整除 < > <= >= == != 小于 大于 小于或等于 大于或等于 等于 不等于 逻辑运算 !x x&y x && y x|y x || y xor(x, y) 逻辑非 逻辑与 同上 逻辑或 同上 异或 数 学 运 算 符 和 比 较 运 算 符 作 用于 两 个 元 素 上(x + y, a < b); 数 学 运 算 符不只是作用于数值型或复数型变量,也可以作用在逻辑型变量上;在后一 种情况中,逻辑型变量被强制转换为数值型。比较运算符可以适用于任何类 型:结果是返回一个或几个逻辑型变量。 逻 辑 型 运 算 符 适 用于 一 个(对 “!” 运 算 符)或 两 个 逻 辑 型 对 象 ( 其 它 运 算 符 ), 并 且 返 回 一 个 ( 或 几 个 ) 逻 辑 性 变 量 。 运 算 符 “ 逻 辑 与 ” 和 “ 逻 辑或”存在两种形式:“&”和“|”作用在对象中的每一个元素上并且返回 和比较次数相等长度的逻辑值;“&&”和“||”只作用在对象的第一个元素 上。 对于0 < x < 1这种类型的不等式必须使用算子“逻辑与”,这个不等式 将 被 改 写 成0 <x & x < 1。 表 达 式0 < x < 1是 合 法 的 , 但 是 不 会 返 回 所 预 期 的 结 果 : 因 为 两 个 运 算 符 相 同 , 它 们 将 从 左 至 右 连 续 执 行 。 首 先 完 成0 < x并 返 回 一 个 逻 辑 变 量 , 它 将 与1比 较(TRUE或 者FALSE< 1): 在这 种 情 形 下 , 这个逻辑变量被强制转换为数值(1或0 < 1)。 > x <- 0.5 >0<x<1 [1] FALSE 10 这些也是R中的运算符:$, @, [, [[, :, ?, <-, <<-, =, ::。可以用?Syntax命令得到描述运 算符之间的优先顺序表。 26 比较运算符作用在两个被比较对象的每个元素上(如果需要,将循环使 用最短的变量),从而返回一个同样大小的对象。为了“整体”比较两个对 象,可以使用两个函数:identical和all.equal > x <- 1:3; y <- 1:3 > x == y [1] TRUE TRUE TRUE > identical(x, y) [1] TRUE > all.equal(x, y) [1] TRUE identical比 较 数 据 的 内 在 关 系 , 如 果 对 象 是 严 格 相 同 的 返 回TRUE, 否 则 返 回FALSE。all.equal用 来 判 断 两 个 对 象 是 否 “ 近 似 相 等 ” , 返 回 结 果 为TRUE或者对二者差异的描述。后一个函数在比较数值型变量时考虑到了计 算过程中的近似。在计算机中数值型变量的比较有时很是令人惊奇。 > 0.9 == (1 - 0.1) [1] TRUE > identical(0.9, 1 - 0.1) [1] TRUE > all.equal(0.9, 1 - 0.1) [1] TRUE > 0.9 == (1.1 - 0.2) [1] FALSE > identical(0.9, 1.1 - 0.2) [1] FALSE > all.equal(0.9, 1.1 - 0.2) [1] TRUE > all.equal(0.9, 1.1 - 0.2, tolerance = 1e-16) [1] "Mean relative difference: 1.233581e-16" 3.5.4 访问一个对象的数值:下标系统 下标系统可以用来有效、灵活且有选择性地访问一个对象中的元素;下 标可以是数值型的或逻辑型的。举例来说,要访问向量x中的第3个值,我们 只要输入x[3]就可以获取或改变这个值: > x <- 1:5 > x[3] [1] 3 > x[3] <- 20 27 >x [1] 1 2 20 4 5 下标本身也可以是一个数值型向量: > i <- c(1, 3) > x[i] [1] 1 20 如果x是一个矩阵或者数据框,第i 行第j 列的值可以通过x[i,j]来访问。 要 访 问 一 个给 定的 行 或 列 中 所 有 的 值 , 只 要 忽 略 适 当的 下 标 ( 不 要 忘 记 逗 号!): > x <- matrix(1:6, 2, 3) >x [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 > x[, 3] <- 21:22 >x [,1] [,2] [,3] [1,] 1 3 21 [2,] 2 4 22 > x[, 3] [1] 21 22 读者必然注意到了最后一个结果是一个向量而不是一个矩阵。R的缺省规 则是返回一个维数尽可能低的对象。这可以通过修改选项drop的值来改变, 它的缺省值是TRUE: > x[, 3, drop = FALSE] [,1] [1,] 21 [2,] 22 下标系统也普遍适用于数组,使用和数组维数同样多的下标(例如,一 个 三 维 数 组 :x[i, j, k], x[, , 3], x[, , 3, drop = FALSE], 等等 ) 。 记住下标必须使用方括号,而圆括号是用来指定函数的参数的: > x(1) Error: couldn’t find function "x" 28 通 过 使 用 负 数 下 标 也 可 以 用 来 不 显 示 一 个 或 一 些行 或 列 。 比 如 ,x[-1,将 不 显 示 第 一 行 ,x[-c(1,15), ]将 不 显 示 第1到第15行 。 使 用 上 面 定 义 的 矩阵: > x[, -1] [,1] [,2] [1,] 3 21 [2,] 4 22 > x[, -(1:2)] [1] 21 22 > x[, -(1:2), drop = FALSE] [,1] [1,] 21 [2,] 22 对于向量,矩阵和数组,可以用一个比较运算表达式作为下标来访问元 素的值: > x <- 1:10 > x[x >= 5] >x [1] 1 2 > x[x == 1] >x [1] 25 2 <- 20 3 4 20 20 20 20 20 20 <- 25 3 4 20 20 20 20 20 20 下面是使用逻辑型下标的一种应用(选择一个整数变量中的偶数): > x <- rpois(40, lambda=5) >x [1] 5 9 4 7 7 6 4 5 11 3 [21] 4 6 6 5 4 5 3 4 3 3 > x[x %% 2 == 0] [1] 4 6 4 2 2 2 4 6 6 4 4 8 4 2 4 5 3 7 7 1 7 5 3 3 8 9 1 2 4 2 2 5 1 2 4 在上面的例子中,下标系统使用了比较运算符返回的逻辑值。这些逻辑 值可以被预先计算,如果需要,将循环使用: > x <- 1:40 > s <- c(FALSE, TRUE) > x[s] [1] 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 29 逻辑下标也可以在数据框中使用,但是要注意数据框中不同的列可能有 不同的数据类型。 对于列表,访问不同的元素(可以是任何类型的对象)可以通过单一的 或者双重的方括号来实现;它们的区别是:单个括号返回一个列表,而双 重括号将提取列表中的对象。在得到的结果中也可以使用下标,就像之前在 向 量 、 矩 阵 等 情 况 中 看 到的 那 样 。 例 如 , 一 个 列 表 中 的第3个 对 象 是 一 个 向 量 , 它 的 取 值 可 以 使 用my.list[[3]][i]来 访 问 , 如 果 是 一 个 三 维 数 组 则 使 用my.list[[3]][i, j, k]等等 ; 另 一 个 区 别 是my.list[1:2]将 返 回 一 个 列 表,包含原始列表的第1个和第2个元素,而my.list[[1:2]]不会给出期望的 结果。 3.5.5 访问对象的名称 names是一个对象元素的字符型标签,它们一般情况下是可选的属性,名 称有多个种类(names, colnames, rownames, dimnames )。 names 是一个和对象有同样长度的向量并且可以通过函数names来访问。 > x <- 1:3 > names(x) NULL > names(x) <- c("a", "b", "c") >x abc 123 > names(x) [1] "a" "b" "c" > names(x) <- NULL >x [1] 1 2 3 对于矩阵和数据框,colnames 和rownames 分别是列和行的标签。它们可 以 通 过 各 自 的 函 数 来 访 问 , 或 者 通 过dimnames 返 回 包 含 两 个 名 称 向 量 的 列 表。 > > > > X <- matrix(1:4, 2) rownames(X) <- c("a", "b") colnames(X) <- c("c", "d") X cd a13 b24 > dimnames(X) 30 [[1]] [1] "a" "b" [[2]] [1] "c" "d" 对于数组,可以用dimnames来访问各维的名字: > A <- array(1:8, dim = c(2, 2, 2)) >A ,,1 [,1] [,2] [1,] 1 3 [2,] 2 4 ,,2 [,1] [,2] [1,] 5 7 [2,] 6 8 > dimnames(A) <- list(c("a", "b"), c("c", "d"), c("e", "f")) >A ,,e cd a13 b24 ,,f cd a57 b68 如果一个对象的元素有名称,它们可以通过名称被提取,和使用下标 一样 。 实 际 上 , 因 为 原 始 对 象 的 属 性 将 仍 被保 留 , 所 以 这 里 的 “ 提 取 ” 应 该 叫 作 “ 取 子 集 ” 更 合 适 。 例 如 , 如 果 数 据 框DF包 含 变 量x,y和z, 命 令DF["x"]将 返 回 一 个 只 包 含x的 数 据 框 ;DF[c("x", "y")]将 返 回 包 含 两 个 变量的数据框。如果列表中的对象有名称,这种方法也将在列表中起作用。 31 就像读者必然认识到的那样,在这里使用的下标是一个字符型向量。就 像在前面看到的数值型或逻辑型向量那样,这个向量可以被预先定义然后用 来提取相应的值。 要 从 一 个 数 据 框 中 提 取 一 个 向 量 或 者 一 个 因 子 , 可 以 使 用 运 算 符$(例 如DF$x),这对列表同样有效。 3.5.6 数据编辑器 我们可以使用一个类似于电子表格的图形编辑器去编辑一个“数据”对 象。例如,如果X是一个矩阵,命令data.entry(X)将打开一个图形编辑器并 且可以通过点击适当的单元格修改数值或者添加新的行或列。 函数data.entry可以直接修改作为其参数给出的对象,而不需要对其结 果赋值。另一方面,函数de可以返回一个列表,其中包含作为函数参数给出 的原始对象(可能已被修改过)。在缺省情况下这个结果被输出在屏幕上, 但是,像对大多数函数一样,它也可以被赋值给一个对象。 使用数据编辑器的具体情况取决于操作系统。 3.5.7 数学运算和一些简单的函数 在R中有很多用来处理数据的函数,前面我们已经看到了最简单的一个函 数c,它用来连接列在圆括号中的对象。如: > c(1:5, seq(10, 11, 0.2)) [1] 1.0 2.0 3.0 4.0 5.0 10.0 10.2 10.4 10.6 10.8 11.0 向量可以进行那些常规的算术运算: >x >y >z >z [1] <- 1:4 <- rep(1, 4) <- x + y 2345 不 同 长 度的 向 量 可 以 相 加 , 这 种 情 况 下 最 短的 向 量 将 被 循 环 使 用 。 例 如: >x >y >z >z [1] >x >y <- 1:4 <- 1:2 <- x + y 2446 <- 1:3 <- 1:2 32 > z <- x + y Warning message: longer object length is not a multiple of shorter object length in: x + y >z [1] 2 4 4 注意R返回了一个警告消息而不是一个错误消息,因此这个操作实际上是 被执行了的。如果我们想要对于一个向量中所有的元素加(或乘)相同的数 值: >x >a >z >z [1] <- 1:4 <- 10 <- a * x 10 20 30 40 R中 用 来 处 理 数 据 的 函 数 太 多 了 而 不 能 全 部 列 在这 里 。 读 者 可 以 找 到 所 有 的 基 本 数 学 函 数(log, exp, log10, log2, sin, cos, tan, asin, acos, atan, abs, sqrt, . . . ), 专业函数(gamma, digamma, beta, besselI, . . . ),同样包括各 种统计学中有用的函数。下表中列出了一部分函数: sum(x) prod(x) max(x) min(x) which.max(x) which.min(x) range(x) length(x) mean(x) median(x) var(x) or cov(x) cor(x) var(x, y) or cov(x, y) cor(x, y) 对x中的元素求和 对x中的元素求连乘积 x中元素的最大值 x中元素的最小值 返回x中最大元素的下标 返回x中最小元素的下标 与c(min(x), max(x))作用相同 x中元素的数目 x中元素的均值 x中元素的中位数 x中元素的的方差(用n − 1做分母);如果x是一个矩阵或者一 个数据框,将计算协方差阵 如 果x是 一 个 矩 阵 或 者 一 个 数 据 框 则 计 算 相 关 系 数 矩 阵 ( 如 果x是一个向量则结果是1) x和y的 协 方 差 , 如 果 是 矩 阵 或 数 据 框 则 计 算x和y对 应 列 的 协 方 差 x和y的 线 性 相 关 系 数 , 如 果 是 矩 阵 或 者 数 据 框 则 计 算 相 关 系 数 矩阵。 这 些 函 数 一 般 返 回 标 量(即 为 长 度 为1的 向 量), 只 有range返 回 一 个 长 度 33 为2的向量,var, cov, 和cor可能会返回一个矩阵。下面的函数返回更复杂的 结果。 round(x, n) rev(x) sort(x) rank(x) log(x, base) scale(x) pmin(x,y,...) pmax(x,y,...) cumsum(x) cumprod(x) cummin(x) cummax(x) match(x, y) which(x == a) choose(n, k) na.omit(x) na.fail(x) unique(x) table(x) table(x, y) subset(x, ...) sample(x, size) 3.5.8 将x中的元素四舍五入到小数点后n位 对x中的元素取逆序 将x中的元素按升序排列;要按降序排列可以用命令rev(sort(x)) 返回x中元素的秩 计算以base为底的x的对数值 如 果x是 一 个 矩 阵 , 则 中 心 化 和 标 准 化 数 据 ; 若 只 进 行 中 心 化 则 使 用 选 项scale=FALSE, 只 进 行 标 准 化 则center=FALSE( 缺 省 值 是center=TRUE, scale=TRUE) 返回一个向量,它的第i个元素是x[i], y[i], . . . 中最小值 同上,取最大值 返回一个向量,它的第i个元素是从x[1]到x[i]的和 同上,取乘积 同上,取最小值 同上,取最大值 返 回 一 个 和x的 长 度 相 同 的 向 量 , 表 示x中 与y中 元 素 相 同 的 元 素 在y中 的位置(没有则返回NA) 返 回 一 个 包 含x符 合 条 件 ( 当 比 较 运 算 结 果 为 真 (TRUE) 的 下 标 的 向 量 , 在这 个 结 果 向 量 中 数 值i说 明x[i] == a( 这 个 函 数 的 参 数 必 须 是 逻辑型变量) 计算从n个样本中选取k个的组合数 忽略有缺失值(NA)的观察数据(如果x是矩阵或数据框则忽略相应的 行) 如果x包含至少一个NA则返回一个错误消息 如果x是一个向量或者数据框,则返回一个类似的对象但是去掉所有重 复的元素(对于重复的元素只取一个) 返回一个表格,给出x中重复元素的个数列表(尤其对于整数型或者因子 型变量) x与y的列联表 返 回x中 的 一 个 满 足 特 定 条 件...的 子 集 , 该 条 件 通 常 是 进 行 比 较 运 算 :x$V1 < 10; 如 果x是数 据 框 , 选 项select给 出 要 保 留 的 变 量 ( 或 者用负号表示去掉) 从x中无放回抽取size个样本,选项replace = TRUE表示有放回的抽样 矩阵计算 R中 有 矩 阵 计 算 和 处 理 的 工 具 , 函 数rbind()和cbind()分 别 用 上 下 或 左 右的方式合并向量或矩阵: > m1 <- matrix(1, nr = 2, nc = 2) 34 > m2 <- matrix(2, nr = 2, nc = 2) > rbind(m1, m2) [,1] [,2] [1,] 1 1 [2,] 1 1 [3,] 2 2 [4,] 2 2 > cbind(m1, m2) [,1] [,2] [,3] [,4] [1,] 1 1 2 2 [2,] 1 1 2 2 两矩阵乘积的运算符是“%*%”,例如,对上面的矩阵m1和m2: > rbind(m1, m2) %*% cbind(m1, m2) [,1] [,2] [,3] [,4] [1,] 2 2 4 4 [2,] 2 2 4 4 [3,] 4 4 8 8 [4,] 4 4 8 8 > cbind(m1, m2) %*% rbind(m1, m2) [,1] [,2] [1,] 10 10 [2,] 10 10 矩阵的转置由函数t完成;这个函数也可以作用于数据框。 函数diag可以用来提取或修正一个矩阵的对角元,或者创建一个对角矩 阵。 > diag(m1) [1] 1 1 > diag(rbind(m1, m2) %*% cbind(m1, m2)) [1] 2 2 8 8 > diag(m1) <- 10 > m1 [,1] [,2] [1,] 10 1 [2,] 1 10 > diag(3) [,1] [,2] [,3] [1,] 1 0 0 35 [2,] 0 1 0 [3,] 0 0 1 > v <- c(10, 20, 30) > diag(v) [,1] [,2] [,3] [1,] 10 0 0 [2,] 0 20 0 [3,] 0 0 30 > diag(2.1, nr = 3, nc = 5) [,1] [,2] [,3] [,4] [,5] [1,] 2.1 0.0 0.0 0 0 [2,] 0.0 2.1 0.0 0 0 [3,] 0.0 0.0 2.1 0 0 R中也有一些用于矩阵计算的专门的函数。我们可以使用solve求矩阵的 逆,用qr来分解矩阵,用eigen来计算特征值和特征向量,用svd来做奇异值 分解。 36 4 R绘图 R提供非常多样的绘图功能。如想了解,可以输入:demo (graphics) 或 者demo(persp)。我们在这里不可能详细说明R在绘图方面的所有功能,主要 是因为每个绘图函数都有大量的选项使得图形的绘制十分的灵活多变。 绘图函数的工作方式与本文前面描述的工作方式大为不同,不能把绘图 函 数 的 结 果 赋给 一 个 对 象11 , 其 结 果 直 接 输 出 到 一 个 “ 绘 图 设 备 ” 上 。 绘 图 设备是一个绘图的窗口或是一个文件。 有 两 种 绘 图 函 数 : 高 级 绘 图 函 数 (high-level plotting functions ) 创 建 一 个新的图形,低级绘图函数(low-level plotting functions )在现存的图形上添 加 元 素 。 绘 图 参 数 (graphical parameters ) 控 制 绘 图 选 项 , 可 以 使 用 缺 省 值 或者用函数par修改。 我们先看如何管理绘图和绘图设备,然后详细说明绘图函数和参数,我 们会看到一个用这些功能产生图形的实例,最后介绍grid和lattice两个包,这 两个包的作法与其他绘图函数不同。 4.1 4.1.1 管理绘图 打开多个绘图设备 当 绘 图 函 数 开 始 执 行 , 如 果 没 有 打 开 绘 图 设 备 , 那 么R将 打 开 一 个 绘 图 窗口来展示这个图形。绘图设备可以用适当的函数打开。可用的绘图设备种 类 取 决 于 操 作 系 统 , 在Unix/Linux 下 , 绘 图 窗 口 称 为x11, 而 在Windows下 称 为windows。 在 所 有 情 况 下 , 都 可 以 用 命 令x11()来 打 开 一 个 绘 图 窗 口 , 在Windows下 仍然 有 效 是 因 为 上 面 的 命 令 可 以 作 为windows()的 别 名 。 可 以 用 函 数 打 开 一 个 文 件 作 为 绘 图 设 备 , 这 包 括 :postscript(), pdf(), png(), . . . 可用的绘图设备列表可以用?device来察看。 最后打开的设备将成为当前的绘图设备,随后的所有图形都在这上面显 示。函数dev.list() 显示打开的列表。 > x11(); x11(); pdf() > dev.list() X11 X11 pdf 2 3 4 显示的数字是设备的编号,要改变当前设备必须使用这些编号,为了解 当前设备用: 11 有一些值得注意的例外: hist()和barplot()仍然生成数据结果作为列表或是矩阵。 37 > dev.cur() pdf 4 为改变当前的设备: > dev.set(3) X11 3 函数dev.off()关闭一个设备:默认关闭当前设备,否则关闭有自变量指 定编号的设备。R然后显示新的当前设备编号。 > dev.off(2) X11 3 > dev.off() pdf 4 在R的Widows版 本 中 , 有 两 个 特 殊 的 功 能 值 得 提 及 :Windows Metafile 设备可以用函数win.metafile来打开,选定绘图窗口会出现“History”菜单, 我 们 可 以 利 用 这 个 菜 单 中 的 功 能 记 录 一 个 会话 中 所 作 的 所 有 图 形 ( 在 缺 省 状 态 下 , 记 录 系 统 是 关 闭 的 , 用 户 可 以 点 击 这 个 菜 单 下 的“Recording”打 开 它)。 4.1.2 图形的分割 函数split.screen分割当前的绘图设备,例如: > split.screen(c(1, 2)) 划分设备为两部分,可以用screen(1)或者screen(2)选择;erase.screen() 删除最后绘制的图形。设备的一部分也可以被split.screen() 划分,可以作 出复杂的布局。 这些函数和其他的函数是不兼容的(比如layout()或者coplot()),不 可以用于多个绘图设备。它们的使用应局限于象图形式探索性数据分析这样 的问题。 函 数layout把 当 前 的 图 形 窗 口 分 割 为 多 个 部 份 , 图 形 将 一 次 显 示 在 各 部分中。它主要的自变量是一个元素都是整数值的矩阵,元素指示子窗口 (“sub-windows”)的编号。例如,把设备划分为4个相等的部分: > layout(matrix(1:4, 2, 2)) 38 当然也可以先产生这个矩阵,以更好的显现设备是如何划分的: > mat <- matrix(1:4, 2, 2) > mat [,1] [,2] [1,] 1 3 [2,] 2 4 > layout(mat) 为了看到创建的分割,我们可以使用函数layout.show,其自变量是子窗 口的个数(这里是4)。在这个例子中,我们有: 1 3 2 4 > layout.show(4) 在下面的例子里,我们将看到layout()提供的各种可能性: 1 5 3 > layout(matrix(1:6, 2, 3)) > layout.show(6) 4 2 > layout(matrix(1:6, 3, 2)) > layout.show(6) 6 1 5 2 > m <- matrix(c(1:3, 3), 2, 2) > layout(m) > layout.show(3) 3 4 6 1 3 2 在 以 上 各个 例 子 中 , 我 们没 有用matrix()的 选 项byrow, 子 窗 口 按 列 编 号 ; 我 们 可 以 指 定matrix(..., byrow=TRUE), 则 窗 口 将 按 行 编 号 。 在 矩 阵 中的编号可以用任何次序,例如matrix(c(2, 1, 4, 3), 2,2)。 缺省情况下,layout()用等间隔分配子窗口:可以用选项widths 和 heights修 改 分 割 的 宽 和 高 。 这 些 尺 寸 是 相 对 给 定 的(也 可 以 用 厘 米 , 详 见?layout),例如: 39 > m <- matrix(1:4, 2, 2) > layout(m, widths=c(1, 3), heights=c(3, 1)) > layout.show(4) 1 3 2 4 2 > m <- matrix(c(1,1,2,1),2,2) > layout(m, widths=c(2, 1), heights=c(1, 2)) > layout.show(2) 1 最后,矩阵里面的编号可以包括0,使得复杂的(甚至怪异的)分割成为 可能。 2 > m <- matrix(0:3, 2, 2) > layout(m, c(1, 3), c(1, 3)) > layout.show(3) 1 > m <- matrix(scan(), 5, 5) 1: 0 0 3 3 3 1 1 3 3 3 11: 0 0 3 3 3 0 2 2 0 5 21: 4 2 2 0 5 26: Read 25 items > layout(m) > layout.show(5) 3 4 1 2 3 5 40 4.2 绘图函数 下面是R中高级绘图函数的概括: plot(x) plot(x, y) sunflowerplot(x, y) pie(x) boxplot(x) stripchart(x) coplot(x~y | z) interaction.plot (f1, f2, y) matplot(x,y) dotchart(x) fourfoldplot(x) assocplot(x) mosaicplot(x) pairs(x) plot.ts(x) ts.plot(x) hist(x) barplot(x) qqnorm(x) qqplot(x, y) contour(x, y, z) filled.contour (x, y, z) image(x, y, z) persp(x, y, z) stars(x) symbols(x, y, ...) termplot(mod.obj) 以x的元素值为纵坐标、以序号为横坐标绘图 x(在x-轴上)与y(在y -轴上)的二元作图 同上 但是以相似坐标的点作为花朵,其花瓣数目为点的个数 饼图 盒形图(“box-and-whiskers”) 把x的值画在一条线段上,样本量较小时可作为盒形图的替代 关于z的每个数值(或数值区间)绘制x与y的二元图 如 果f1和f2是 因 子 , 作y的 均 值 图 , 以f1的 不 同 值 作 为x轴 , 而f2的 不 同 值 对 应 不 同 曲 线 ; 可 以 用 选 项fun指 定y的 其 他 的 统 计量(缺省计算均值,fun=mean) 二元图,其中x的第一列对应y的第一列,x的第二列对应y的第二 列,依次类推。 如果x是数据框,作Cleveland点图(逐行逐列累加图) 用 四 个 四 分 之 一 圆 显 示2X2列 联 表 情 况 (x必 须 是dim=c(2, 2, k)的数组,或者是dim=c(2, 2)的矩阵,如果k = 1) Cohen–Friendly图,显示在二维列联表中行、列变量偏离独立性 的程度 列联表的对数线性回归残差的马赛克图 如果x是矩阵或是数据框,作x的各列之间的二元图 如 果x是 类"ts"的对 象 , 作x的 时 间 序 列 曲 线 ,x可 以 是 多 元 的 , 但是序列必须有相同的频率和时间 同 上 , 但 如 果x是 多 元 的 , 序 列 可 有 不 同 的 时 间 但 须 有 相 同 的 频 率 x的频率直方图 x的值的条形图 正态分位数-分位数图 y对x的分位数-分位数图 等 高 线 图 ( 画 曲 线 时 用 插 补 充 空 白 的 值 ) ,x和y必 须 为 向 量 ,z必 须 为 矩 阵 , 使 得dim(z)=c(length(x), length(y))(x和y可以省略) 同上,等高线之间的区域是彩色的,并且绘制彩色对应的值的图 例 同上,但是实际数据大小用不同色彩表示 同上,但为透视图 如果x是矩阵或者数据框,用星形和线段画出 在由x和y给定坐标画符号(圆,正方形,长方形,星,温度计式 或者盒形图),符号的类型、大小、颜色等由另外的变量指定 回归模型(mod.obj)的(偏)影响图 41 每一个函数,在R里都可以在线查询其选项。某些绘图函数的部分选项是 一样的;下面列出一些主要的共同选项及其缺省值: add=FALSE axes=TRUE type="p" xlim=, ylim= xlab=, ylab= main= sub= 4.3 如果是TRUE,叠加图形到前一个图上(如果有的话) 如果是FALSE,不绘制轴与边框 指定图形的类型,"p": 点,"l": 线,"b": 点连线,"o": 同上,但是线在点上,"h": 垂直线,"s": 阶梯式,垂直 线顶端显示数据,"S": 同上,但是在垂直线底端显示数 据 指 定 轴 的 上 下 限 , 例 如xlim=c(1, 10)或 者xlim=range(x) 坐标轴的标签,必须是字符型值 主标题,必须是字符型值 副标题(用小字体) 低级绘图命令 R里 面 有 一 套 绘 图 函 数 是 作 用 于 现 存 的 图 形 上 的 : 称 为 低 级 作 图 命 令(low-level plotting commands) 。下面有一些主要的: points(x, y) lines(x, y) text(x, y, labels, ...) mtext(text, side=3, line=0, ...) segments(x0, y0, x1, y1) arrows(x0, y0, x1, y1, angle= 30, code=2) abline(a,b) abline(h=y) abline(v=x) abline(lm.obj) rect(x1, y1, x2, y2) polygon(x, y) legend(x, y, legend) title() 添加点(可以使用选项type=) 同上,但是添加线 在(x,y)处添加用labels指定的文字;典型的用法是: plot(x, y, type="n"); text(x, y, names) 在边空添加用text指定的文字,用side指定添加到哪一边(参照 下面的axis());line指定添加的文字距离绘图区域的行数 从(x0,y0)各点到(x1,y1)各点画线段 同 上 但 加 画 箭 头 , 如 果code=2则 在 各(x0,y0)处 画 箭 头 , 如 果code=1则在 各(x1,y1)处 画 箭 头 , 如 果code=3则在 两 端都 画 箭 头; angle控制箭头轴到箭头边的角度 绘制斜率为b和截距为a的直线 在纵坐标y处画水平线 在横坐标x处画垂直线 画由lm.obj确定的回归线(参照第五章) 绘制长方形,(x1, y1)为左下角,(x2,y2)为右上角 绘制连接各x,y坐标确定的点的多边形 在点(x,y)处添加图例,说明内容由legend给定 添加标题,也可添加一个副标题 42 axis(side, vect) box() rug(x) locator(n, type="n", ...) 画坐标轴,side=1时画在下边,side=2时画在左边,side=3时画 在上边,side=4时画在右边。可选参数at指定画刻度线的位置坐 标 在当前的图上加上边框 在x-轴上用短线画出x数据的位置 在 用 户 用 鼠 标 在 图 上 点 击n次 后 返 回n次 点 击 的 坐 标(x, y ); 并 可 以 在 点 击 处 绘 制 符 号(type="p"时)或 连 线(type="l"时), 缺 省 情 况下不画符号或连线 注 意 , 用text(x, y,expression(...))可 以 在 一 个 图 形 上 加 上 数 学 公 式,函数expression把自变量转换为数学公式。例如, > text(x, y, expression(p == over(1, 1+e^-(beta*x+alpha)))) 在图中相应坐标点(x, y )处显示下面的方程: p= 1 1+ e−(β x+α) 为了能在表达式中代入某个变量的值,我们可以使用函数substitute和 as.expression,例如,为了代入R2 的值(之前计算并储存在对象Rsquared中) > text(x, y, as.expression(substitute(R^2==r, list(r=Rsquared)))) 在图中相应坐标点(x, y )处显示: R2 = 0.9856298 如果只显示3位小数,我们可以修改代码如下: > text(x, y, as.expression(substitute(R^2==r, + list(r=round(Rsquared, 3))))) 将显示: R2 = 0.986 最后,用斜体字显示R: > text(x, y, as.expression(substitute(italic(R)^2==r, + list(r=round(Rsquared, 3))))) R2 = 0.986 43 4.4 绘图参数 除了低级作图命令之外,图形的显示也可以用绘图参数来改良。绘图参 数可以作为图形函数的选项(但不是所有参数都可以这样用),也可以用函 数par来永久地改变绘图参数,也就是说后来的图形都将按照par指定的参数 来绘制。例如,下面的命令: > par(bg="yellow") 将导致后来的图形都以黄色的背景来绘制。有73个绘图参数,其中一些 有非常相似的功能。这些参数详细的列表可以参阅?par;下面的表格只列举 了最常用的参数。 adj bg bty cex col font las lty lwd mar mfcol mfrow pch ps pty tck tcl xaxt yaxt 控制关于文字的对齐方式,0是左对齐,0.5是居中对齐,1是右对齐,值> 1时 对齐位置在文本右边的地方,取负值时对齐位置在文本左边的地方;如果给出 两个值(例如c(0, 0)),第二个只控制关于文字基线的垂直调整 指 定 背 景 色 ( 例 如bg="red", bg="blue"; 用colors()可 以 显 示657种 可 用 的 颜 色名) 控制图形边框形状,可用的值为: "o", "l", "7", "c", "u" 和"]" (边框和字符 的外表相像);如果bty="n"则不绘制边框 控 制 缺 省 状 态 下 符 号 和 文 字 大 小 的 值 ; 另 外 ,cex.axis控 制 坐 标 轴 刻 度 数 字 大 小 ,cex.lab控 制 坐 标 轴 标 签 文 字 大 小 ,cex.main控 制 标 题 文 字 大 小,cex.sub控制副标题文字大小 控 制 符 号 的 颜 色 ; 和cex类 似 , 还 可 用 :col.axis, col.lab, col.main, col.sub 控制文字字体的整数(1: 正常,2: 斜体,3: 粗体,4: 粗斜体);和cex类似, 还可用: font.axis, font.lab, font.main, font.sub 控制坐标轴刻度数字标记方向的整数(0: 平行于轴,1: 横排,2: 垂直于轴,3: 竖排) 控 制 连 线 的 线 型 , 可 以 是 整 数 (1: 实 线 ,2: 虚 线 ,3: 点 线 ,4: 点 虚 线 ,5: 长 虚 线 ,6: 双 虚 线 ) , 或 者 是 不 超 过8个 字 符 的 字 符 串 ( 字 符 为 从"0"到"9"之 间 的 数 字 ) 交 替 地 指 定 线 和 空 白 的 长 度 , 单 位为 磅(points)或 象素,例如lty="44"和lty=2效果相同 控制连线宽度的数字 控 制 图 形 边 空 的 有4个 值 的 向 量c(bottom, left, top, right), 缺 省 值 为c(5.1, 4.1, 4.1, 2.1) c(nr,nc)的向量,分割绘图窗口为nr行nc列的矩阵布局,按列次序使用各子窗 口(参照 4.1.2) 同上,但是按行次序使用各子窗口(参照 4.1.2) 控制符号的类型,可以是1到25的整数,也可以是""里的单个字符(图 2) 控制文字大小的整数,单位为磅(points) 指定绘图区域类型的字符,"s": 正方形,"m":最大利用 指定轴上刻度长度的值,单位是百分比,以图形宽、高中最小一个作为基数; 如果tck=1则绘制grid 同上,但以文本行高度为基数(缺省下tcl=-0.5) 如果xaxt="n"则设置x-轴但不显示(有助于和axis(side=1, ...)联合使用) 如果yaxt="n"则设置y -轴但不显示(有助于和axis(side=2, ...)联合使用) 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 "*" "?" "." "X" "a" * ? Xa 图 2: R (pch=1:25)的绘图符号。用选项col="blue", bg="yellow"来产生如 上的颜色,其中背景色选项只对符号21–25有作用。可以使用任意字符作为绘 点符号(pch="*", "?", ".", . . . )。 4.5 一个实例 为了讲解R的绘图功能,让我们来看一个简单的10对随机值的二维图形的 例子。这些值用以下命令生成: > x <- rnorm(10) > y <- rnorm(10) 所需的图可以用plot()来产生;只要输入命令: > plot(x, y) 则图形将绘制在当前的绘图设备上。结果见图 3。缺省情况下,R用“智 能”的方法绘制图形:R自动计算坐标轴上刻度摆放,标记的位置等,使得图 形尽可能的易于理解。 用户仍然可以改变绘图的法,例如,为了遵照某刊物的要求,或为某 个演讲作个性化调整。最简单的方式就是用选项值取代缺省值来修改图形绘 制方式。在我们的例子中,我们可以用以下方式大大改变图形: plot(x, y, xlab="Ten random values", ylab="Ten other values", xlim=c(-2, 2), ylim=c(-2, 2), pch=22, col="red", bg="yellow", bty="l", tcl=0.4, main="How to customize a plot with R", las=1, cex=1.5) 45 0.5 0.0 −1.0 −0.5 y −0.5 0.0 0.5 1.0 x 图 3: 没有用任何选项的函数plot。 结果见图 4。我们来详细说明其中的每个选项。首先,xlab 和ylab改变 坐 标 轴 标 签 , 缺 省 情 况 下 是 变 量 的 名 字 。 然 后 ,xlim 和ylim允 许 我 们 规 定 两个坐标轴的范围12 。绘图参数pch在这里用作一个选项:pch=22为正方形, 其轮廓颜色和背景色可能不一样,分别由col 和bg指定。图形参数表中说明 了bty,tcl,las和cex的作用最后,选项main添加了标题。 绘图参数和低级作图函数使我们可以进一步改善图形。前面我们已经看 到,一些绘图参数不允许作为plot这样的函数的自变量。我们可以用par()来 修改这些参数,这样就必须输入多行的命令。在改变绘图参数时,预先保存 它们的初始值以便以后恢复十分有用。以下命令可以产生图 5。 opar <- par() par(bg="lightyellow", col.axis="blue", mar=c(4, 4, 2.5, 0.25)) plot(x, y, xlab="Ten random values", ylab="Ten other values", xlim=c(-2, 2), ylim=c(-2, 2), pch=22, col="red", bg="yellow", bty="l", tcl=-.25, las=1, cex=1.5) title("How to customize a plot with R (bis)", font.main=3, adj=1) par(opar) 我 们 详细 解 释 这 些 命 令 的 作 用 。 首 先 , 缺 省 的 绘 图 参 数 被 复 制 到 列 表opar中。然后有三个参数被修改:bg修改背景色,col.axis修改轴上数字 的颜色,mar来修改绘图边空大小。这个图形用几乎和图 4相似的方式画出。 边 空 的 修 改 让 我 们 更 好 地 利 用 绘 图 区 周 围 的 空 白 。 用 低 级 函 数title添 加 标 12 缺 省 下 ,R在 每 个 坐 标 轴 界 限 上 增 加4%。 这 种 习 惯 可 以 随 着 设 置 绘 图 参 数xaxs="i"和yaxs="i"而改变(他们可以作为plot()的选项被传递)。 46 How to customize a plot with R 2 Ten other values 1 0 −1 −2 −2 −1 0 1 2 Ten random values 图 4: 用了选项的函数plot。 How to customize a plot with R (bis) 2 Ten other values 1 0 −1 −2 −2 −1 0 1 Ten random values 图 5: 函数par,plot和title。 47 2 题,这样允许给定一些参数作为自变量而不改变图形的其余部分。最后一行 的命令恢复初始的绘图参数。 现在,完全控制!在图 5中,R仍然自动决定了诸如坐标轴刻度的个数, 标题与绘图区域之间的距离等少许事情。我们现在将看到如何完全控制图形 的绘制。这里用的方法是用plot(...,type="n")绘制一个“空白”的图形, 然后用低级函数来添加点,坐标轴,标签等。我们可以想出诸如改变绘图区 域颜色这样的安排。命令如下,产生的图形见图 6。 opar <- par() par(bg="lightgray", mar=c(2.5, 1.5, 2.5, 0.25)) plot(x, y, type="n", xlab="", ylab="", xlim=c(-2, 2), ylim=c(-2, 2), xaxt="n", yaxt="n") rect(-3, -3, 3, 3, col="cornsilk") points(x, y, pch=10, col="red", cex=2) axis(side=1, c(-2, 0, 2), tcl=-0.2, labels=FALSE) axis(side=2, -1:1, tcl=-0.2, labels=FALSE) title("How to customize a plot with R (ter)", font.main=4, adj=1, cex.main=1) mtext("Ten random values", side=1, line=1, at=1, cex=0.9, font=3) mtext("Ten other values", line=0.5, at=-1.8, cex=0.9, font=3) mtext(c(-2, 0, 2), side=1, las=1, at=c(-2, 0, 2), line=0.3, col="blue", cex=0.9) mtext(-1:1, side=2, las=1, at=-1:1, line=0.2, col="blue", cex=0.9) par(opar) 和以前一样,先保存缺省的绘图参数,然后修改背景颜色和边空。画图 时用type="n"不画出点,用xlab="",ylab=""不画坐标轴标签,和用xaxt="n" ,yaxt="n"不 画 坐 标 轴 。 这 样 只 画 了 绘 图 区 域 的 边 框 , 并 用xlim和ylim规 定 了坐标轴范围。注意,我们可以用选项axes=FALSE,但这样的话不仅不画坐 标轴,而且也不画边框。 然后,用低级图形函数在上面确定的坐标区域内加入各种图形元素。在 添加点以前,用rect()修改绘图区域的颜色:长方形大小选得比绘图区域大 得多。 用points()画点;用了一个新的符号。用axis()添加坐标轴:第二个自 变量提供的向量指定坐标刻度位置。选项labels=FALSE指定画坐标轴时不画 刻度数字。这个选项也可以用于字符式样的向量,例如labels=c("A", "B", "C")。 用title()添 加 标 题 , 但 是 字 体 稍 微 改 变 了 。 开 始 的 两 个 边 空 文 字 函 数mtext()调 用 画 坐 标 轴 的 标 签 。 这 个 函 数 的第 一 自 变 量 是 要 画 的 文 本 。 选 项line指出到绘图区域的距离行数(缺省时line=0),at给出坐标。第二次 48 Ten other values How to customize a plot with R (ter) 1 0 −1 −2 0 Ten random values 2 图 6: “手工”图。 调用mtext(),调用利用了side (3)的缺省值。另外两个mtext()用数值型向量 作第一自变量,会自动转换为字符型。 4.6 grid 和lattice 包 grid 和lattice包实现grid和lattice系统。grid是一个新的绘图模式,其绘图 参数的系统和上面所讲的截然不同。grid与基本的图形比较,主要的不同之处 有以下两点: • 用viewports可以以更灵活的方式来分割绘图设备,使其全部覆盖(绘图 对象甚至可以在不同的视口中共享,例如,箭头记号); • 绘图对象(grob)可以被修改或者从一个图形中移除,并不需要重新绘制 所有的图形(用基本图形时必须重绘)。 grid图形通常不可以用来和基本图形组合或者混合(必须用gridBase包来 作 ) 。 可 是 , 在 相 同 的 绘 图 设 备 上 , 同 一 会话 里 用 两 种 绘 图 模 式是 有 可 能 的。 Lattice本来是S-PLUS里的Trellis图形在R中的实现。Trellis是多元数据可 视化的方法,特别适用于发现各变量之间的相互作用关系13 。Lattice(Trellis) 的主要想法是不同条件下的多个图:根据某变量的值的不同对两个变量作不 同图。函数coplot作用类似,但是lattice提供了更广泛的功能。Lattice用grid图 形模式。 13 http://cm.bell-labs.com/cm/ms /departments/sia/project/trellis/index.html 49 在lattice里 , 大 部 分 的 函 数 都 把 一 个公 式 作 为 他 们 的 主 要 自 变 量14 , 例 如y ~ x。公式y ~ x | z是指按z值的不同关于y对x作不同图形。 下 面 的 表 格给 出 了lattice中 的 主 要 函 数 。 表 中 列 出 的 作 为 自 变 量 的 公 式 是典型的用法,但是,所有的这些函数接受条件公式(y ~ x | z)作为主要变 量;在后一种情况下,对z的不同取值将作多个图形。 barchart(y ~ x) bwplot(y ~ x) densityplot(~ x) dotplot(y ~ x) histogram(~ x) qqmath(~ x) stripplot(y ~ x) qq(y ~ x) xyplot(y ~ x) levelplot(z ~ x*y) contourplot(z ~ x*y) cloud(z ~ x*y) wireframe(z ~ x*y) splom(~ x) parallel(~ x) y对x的直方图 盒 形图 密度函数图 Cleveland点图(逐行逐列累加图) x的频率直方图 x的关于某理论分布的分位数-分位数图 一维图,x必须是数值型,y可以是因子 比 较 两 个 分 布 的 分 位 数 ,x必 须 是数 值 型 ,y可 以 是数 值 型,字符型,或者因子,但是必须有两个”水平” 二元图(有许多功能) 在x,y坐标点的z值的彩色等值线图(x,y和z等长) 3-D透视图(点) 同上(面) 二维图矩阵 平行坐标图 为 了 举 例 说 明lattice的 一 些 方 面 , 让 我 们 来 看 一 些 例 子 。 包 必 须 要 用 命 令library(lattice)来装载在内存中,使得包的函数可以被访问。 让我们从密度函数图形开始。这样的图形可以用densityplot(~ x)简单 的作出,将作出经验密度函数曲线并在x-轴处用散点显示各观测值(如rug()所 作)。我们的例子将稍微的变复杂一些,在每个图形里,除经验密度曲线之外 还叠加一个正态密度拟合曲线。这样必须用自变量panel 来定义每个图上绘 制什么。命令如下: n <- seq(5, 45, 5) x <- rnorm(sum(n)) y <- factor(rep(n, n), labels=paste("n =", n)) densityplot(~ x | y, panel = function(x, ...) { panel.densityplot(x, col="DarkOliveGreen", ...) panel.mathdensity(dmath=dnorm, args=list(mean=mean(x), sd=sd(x)), col="darkblue") }) 14 plot() 也 接 受 一 个公 式 作 为 其 主 要 变 量 : 如 果x和y是 两 个 相 同 长 度的 向 量 ,plot(y ~ x)和plot(x,y)将绘制同样的图形。 50 −4 −2 0 2 4 n = 35 n = 40 n = 45 n = 20 n = 25 n = 30 0.6 0.5 0.4 0.3 0.2 0.1 0 0.6 Density 0.5 0.4 0.3 0.2 0.1 0 n=5 n = 10 n = 15 0.6 0.5 0.4 0.3 0.2 0.1 0 −4 −2 0 2 4 −4 −2 0 2 4 x 图 7: 函数densityplot。 命 令 的 前 三 行 产 生 随 机 独 立 正 态 样 本 , 分 割 成 个 数 等 于5, 10, 15, . . . , 和45 的子样本。然后用densityplot为每个子样本产生图形。panel作为函数 的自变量。在我们的例子中,我们定义一个函数调用lattice中的预先确定的两 个 函 数 :panel.densityplot绘 制 经 验 密 度 函 数 ,panel.mathdensity绘 制 拟 合 的 正 态 分 布 密 度 函 数 。 函 数panel.densityplot被 缺 省 调 用 如 果 没 有 自 变量给panel:命令densityplot (~ x | y)将有和图 7相同的结果,但是没 有蓝线。 下一个例子是修改lattice包帮助理的例子得到的,用一些R里面有的数据 集:斐济岛附近1000个地震位置,和一些三种鸢尾花朵的测量值。 图 8显示不同深度的地震的地理位置,作图命令如下: data(quakes) mini <- min(quakes$depth) maxi <- max(quakes$depth) int <- ceiling((maxi - mini)/9) inf <- seq(mini, maxi, int) quakes$depth.cat <- factor(floor(((quakes$depth - mini) / int)), labels=paste(inf, inf + int, sep="-")) xyplot(lat ~ long | depth.cat, data = quakes) 第一个命令是在内存中装载数据quakes。后面的5行命令创建了一个把深 度等分为九个区间的因子:因子水平标签为区间的上下界。下面就只要用适 当的公式调用xyplot函数,其中用data自变量指定绘图用的各变量所在的数 51 165 170 175 180 185 472−544 544−616 616−688 256−328 328−400 400−472 −10 −15 −20 −25 −30 −35 −40 −10 −15 lat −20 −25 −30 −35 −40 40−112 112−184 184−256 −10 −15 −20 −25 −30 −35 −40 165 170 175 180 185 165 170 175 180 185 long 图 8: 应用数据“quakes”的函数xyplot。 据框15 。 关 于 数 据iris, 不 同 物 种 的 重 叠 部 分 十 分 小 , 使 得 可 以 画 在 图 9中 。 命 令如下: data(iris) xyplot( Petal.Length ~ Petal.Width, data = iris, groups=Species, panel = panel.superpose, type = c("p", "smooth"), span=.75, auto.key = list(x = 0.15, y = 0.85) ) 在这 里 调 用 函 数xyplot比 之 前 的 例 子 要 复 杂 些 , 所 用 的 一 些选 项 , 我 们 将 详 细 介 绍 。 选 项groups定 义 组 以 便 被 其 他 选 项 引 用 。 我 们 已 经 知 道 选 项panel是 来 定 义 如 何 在 图 中 显 示 不 同 组 的 : 我 们 在这 里 用 了 预 先 定 义 的 函 数panel.superpose, 是 为 了 把 不 用 的 组 重 叠的 绘 制 在 同 一 个 图 里 。 没 有 给panel.superpose传 递 选 项 , 使 用 缺 省 颜 色 用于 区 分 不 同 的 组 。 选 项type和在plot()中一样,指定数据如何显示,但是这里我们可以给出包含 若干元素的向量作为自变量:"p" 表示画点,"smooth"表示画光滑曲线,其 光滑程度用span指定。选项auto.key在图形中添加图例:只需要给出一个指 定图例的坐标的列表。注意,这里的坐标是相对与图形大小的比例值(也就 是在[0, 1])。 15 plot() 不 可 以 用data作 为 自 变 量,变 量 的 位 置 必 须 明 确 的 给 定 , 例 如plot(quakes$long ~quakes$lat)。 52 7 o o o setosa versicolor virginica 6 o o o o o o 5 Petal.Length o o o 4 o o o 3 o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o 2 o o o 1 0 o o o o o o o o o o o o o o o o o o 0.5 o 1 1.5 2 2.5 Petal.Width 图 9: 应用数据“iris”的函数xyplot。 我 们 现 在 看 到 用 函 数splom用于 相 同 的 数 据iris上 。 用 以 下 命 令 来 绘 制 图 10: splom( ~iris[1:4], groups = Species, data = iris, xlab = "", panel = panel.superpose, auto.key = list(columns = 3) ) 这 次 主 要 的 自 变 量 是 一 个 矩 阵 (iris的 前 四 列 ) 。 结 果 是 一 系 列 可 能 的 在 矩 阵 列 之中 的 二 元 图 , 就 像 标 准 函 数pairs。 缺 省 下 ,splom在x-下 加 上“Scatter Plot Matrix”的文字:要消除的话,可以用选项xlab=""。其他的 选项和前面的例子类似,除了auto.key的columns = 3是指定图例显示在3列 中。 图 10可 以 用pairs()做 出 , 但 是 后 者 的 函 数 不 可 以 作 像 图Fig. 11一样 的 有条件的图形。使用的代码相对简单: splom(~iris[1:3] | Species, data = iris, pscales = 0, varnames = c("Sepal\nLength", "Sepal\nWidth", "Petal\nLength")) 子图相对比较小,我们添加了两个选项来改进图的易理解性:pscales = 0 移除坐标轴的刻度(所有的子图用相同的刻度绘制),变量的名字被重新 定义,并用两行来显示("\n"表示换行)。 最 后 的 例 子 用 平 行 坐 标 来 探 索 多 元 数 据 分 析 的 方法 。 变 量 安 排 在 一 个 坐 标 轴 上 ( 例 如y -轴 ), 观 测 值 绘 制 在 另 一 坐 标 轴 上 ( 变 量 取 值 范 围 调 整 到 53 Setosa Versicolor oo oo o oo oo o o o oo oo o o o ooo o o oo o o o o o o oo o o oo o o o oo o oo oooo o oo o oo o oo o ooo o o o oo o o o oooooooo oo o o o ooo o oooooooo oo o o o o oo o oo o oo o oo o oooo oo o o o oooooo o oo o o o o o oo oo o o oo o oooooo Virginica oo oooo o o oo o oo o oo ooo o oo o oo o o oo o o o oo o oo ooo oo o o oo oo o oo o ooo o ooo ooooo o ooo o oo oo oo o ooo o 2.5 1.5 o oo o ooo oo o ooo oo ooo oo oo o o oooooooooo ooo oooo oo 0 oo o o o ooo o o 7 o oo 4567 o o o o o oo ooo o o o oo o o o 6 oo oo ooooooo o o o oo o o oo oo o o oo ooo o ooo o 5 o o oo ooo o oo o oo o o oo o o o o oo o o oo oo oo ooooooo o oo o o oo oooo o o o oo o oo oo oo ooooooo o 4 Petal.Length 4 o o o oo oo oo o oo o o o o o o o 3 7 Sepal.Length 6 5 5 6 oo o oooo oo oooooooooooooo o 1 o o o 4.5 3.5 4 3.5 Sepal.Width 2.5 2 2.5 2 2 4.5 o o oo ooo o o o o ooo o o oo o oo o oo oo o o oo oo 3 oo o 4 3 o 2 o o ooo o oo oo oo o ooo o o oooooo ooooo ooo o o o oo o oooo oo o ooooo o o oo oo o o o oo o oooooo ooo oooo oo o o o o oo o o oo o o o o o o ooo o o o ooooooo o o o ooo oo o o oo o oo oo o o o ooo o o o oo oo o o o o ooo o o oo o o o o o 3 4 1 2 2.5 Petal.Width 1 o o o oo o o o o o oo o oo o oo o o oo o o o oo o o o oo o o o oo o oo o oo o oo o o o oo o o o o o o ooo ooo oo o o o oo o oo o o o o o oo o o o o o o o oo o o ooo o o oo o o o oo o o o o o o oo o o o o o o o o oo o oo o o o o o o oo o oo o o ooo o o o o oo o ooo o o oo oo o oo oo oo oo o o o o o oo oo o o 8 7 8 1.5 2 0.5 0.5 o oo oooooo oooo oo o o o oo oo o ooo oooo o o oooo o o ooo oo o oo o oo o oo oo o oo o oo o oo oo oooo o oo oooo o o o o oooooooo o o oooo o o o ooo oo o o ooo o o o o oo oo oo oo o o o o oo oo o o o o o o o oo o oo oo o o ooo oo o o ooo o o oo ooo oo o oo oo ooo oo o o oo o ooo o o o oo o o ooo oo o ooo o oo o ooo o o oooo o oo oo o oo o o o o ooooo oo o o o oo oooo o ooo o oo oo o 1 0 oooo ooo o o ooo o o o ooooo o oo oo oo oooooo o o ooooooo ooooo oo oooooo o oo o o oo o o oo ooo o oo o o oo o ooooooo ooooo ooo o o o ooo oo o ooooo ooooooo o ooo oo o o oooo oo o o o o oooo o oo oo o o oo oo o oo ooo oo oooooo oo o oo o oo oo ooooo o o oo o o o o ooooo o o oo o ooo o oooo o o o oooo o oo o o o 图 10: 应用数据数据“iris”(1)的函数splom。 virginica Petal Length Sepal Width Sepal Length setosa versicolor Petal Length Petal Length Sepal Width Sepal Width Sepal Length Sepal Length Scatter Plot Matrix 图 11: 应用数据“iris”(2)的函数splom。 54 Min setosa Max versicolor virginica Petal.Width Petal.Length Sepal.Width Sepal.Length Min Max Min Max 图 12: 应用数据“iris”的函数parallel。 可比,例如,标准化之)。相同个体不同变量的值连接在一条线上。用数 据iris,下面的命令可以得到图 12: parallel(~iris[, 1:4] | Species, data = iris, layout = c(3, 1)) 55 5 R的统计 分析 尽管R的统计功能比作图功能强很多,但这里不可能很详细地去描述它所 有的统计分析功能. 我的目的是在这里给读者对R统计分析功能的一个粗略而 又系统的介绍. 包stats 包 括 了 一 系 列 基 本 的 统 计 分 析 函 数 : 经 典的 假 设 检 验, 线 性 模 型(包括最小二乘法回归, 广义线性模型, 和方差分析), 统计分布, 汇总统计, 层 次聚类16 , 时间序列分析, 非线性最小二乘法和多元分析. 其他的R包还提供了 一些上述统计方法以外的的统计方法. 和基本R安装同时发布的统计包被标注 为推荐包, 而其他的包标注为捐献包并且要求用户自己安装. 我 们 以一 个 简 单 的 例 子 开 始 。 这 个 例 子 仅仅 需 要 包stats, 但 它 可 以 说 明 用R进行分析的大体过程。然后,我们将细致讲述两个在所有统计分析中都非 常有用的概念, 公式(formulae) 和泛型函数(generic functions). 我们还会对所 有的包进行一个总结. 5.1 关于方差分析的一个简单例子 在 包stats 里 面 的 方 差 分 析 函 数是aov. 为 了 示 范 这 个 函 数, 我 们 采 用R分 发的数据集: InsectSprays. 六种杀虫剂将会在田野中进行效价测试, 观测变 量是昆虫的个数. 每种杀虫剂重复测试了12 次, 因此共有72次观测致.我们不 想在这里讨论这些数据的统计图分析, 我们的焦点在于研究响应变量关于杀虫 剂种类的简单方差分析. 通过函数data 把数据导入内存后, 首先对数据进行平 方根转换,然后再进行分析: > data(InsectSprays) > aov.spray <- aov(sqrt(count) ~ spray, data = InsectSprays) 函数aov的主要参数(必要的)是一个公式. 公式的左边是响应变量,右边 是 预 测变 量 , 二 者 通 过~ 连 接. 可 选 项data = InsectSprays 表 明 这 些 变 量 是在数据框InsectSprays中. 下面的语句和前面的语句在语法上等价: > aov.spray <- aov(sqrt(InsectSprays$count) ~ InsectSprays$spray) 如果我们知道列的编号, 我们还可以如下分析: > aov.spray <- aov(sqrt(InsectSprays[, 1]) ~ InsectSprays[, 2]) 16 译者注:hierarchical clustering 56 当然, 第一种方式是首选的, 因为它最为简明易懂. 前面脚本运行的结果不会在屏幕上显示, 因为它们都被赋给对象aov.spray. 我们必须采用一些函数去解析这些结果, 如函数print 可以对分析结果进行一 个简单的总结(一般是要估计的参数), 函数summary 可以显示更多的细节(包括 统计检验的结果): > aov.spray Call: aov(formula = sqrt(count) ~ spray, data = InsectSprays) Terms: spray Residuals Sum of Squares 88.43787 26.05798 Deg. of Freedom 5 66 Residual standard error: 0.6283453 Estimated effects may be unbalanced > summary(aov.spray) Df Sum Sq Mean Sq F value Pr(>F) spray 5 88.438 17.688 44.799 < 2.2e-16 *** Residuals 66 26.058 0.395 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 需要说明的是,直接键入对象名字和命令print(aov.spray) 的结果是类 似的. 可以通过函数plot() 和termplot() 展示分析结果的统计图. 在键入函 数plot(aov.spray)前, 我们必须把图界面分成四部分, 这样就可以让四个 统计诊断图在同一个图上显示. 相应的命令是17 : > > > > > opar <- par() par(mfcol = c(2, 2)) plot(aov.spray) par(opar) termplot(aov.spray, se=TRUE, partial.resid=TRUE, rug=TRUE) 运行结果如图 13 和14 所示. 17 译者注:我个人更喜欢用函数layout 57 1.5 2.5 0.5 3.5 1.5 1 2 0.08 27 25 0 Theoretical Quantiles 39 0.04 Cook’s distance Cook’s distance plot 25 0 3.5 0.00 0 1 2 39 27 −1 2.5 Fitted values Normal Q−Q plot −2 Standardized residuals Fitted values −2 39 1.0 1.5 Scale−Location plot 27 25 0.0 25 Standardized residuals 1.0 39 0.0 −1.5 Residuals Residuals vs Fitted 27 20 40 60 Obs. number 图 13: 通过函数plot() 图形化展示函数aov的分析结果. 5.2 公式 公式是R统计分析里面的关键元素: 几乎所有函数都采用一样的符号18 .公 式的典型形式是y ~ model , 其中y是响应变量, model 是一些元素项的集合而 且要为其中一些项估计参数. 这些元素项通过一些有特殊涵义的运算符连接. a+b X a:b a*b poly(a, n) ^n b %in% a -b -1 1 offset(...) 18 a 和b的相加效应 如 果X 是 一 个 矩 阵, 这 将 反 映 各 列 的 相 加 效 应, 即 X[,1]+X[,2]+...+X[,ncol(X)]; 还 可 以 通 过 索 引 向 量 选择特定列进行分析(如, X[,2:4]) a 和b 的交互效应 相加和交互效应(等价于a+b+a:b) a的n价多项式 包 含 所 有 的 直 到n阶 的 交 互 作 用, 即 (a+b+c)^2 等 价 于a+b+c+a:b+a:c+b:c b 和a的嵌套分类设计(等价于a+a:b, 或者a/b) 去 掉 因 子b的 影 响, 如: (a+b+c)^2-a:b 等 价 于a+b+c+a:c+b:c y~x-1 表 示 通 过 原 点 的 线 性 回 归(等 价 于 y~x+0 或 者0+y~x) y~1 拟合一个没有因子影响的模型(仅仅是截距) 向 模 型 中 增 加 一 个 影 响 因 子 但 不 估 计 任 何 参 数(如, offset(3*x)) 译者注:也有些例外, 基本上趋同. 58 2 1 0 −1 −3 −2 Partial for spray 0 1 2 3 4 5 6 spray 图 14: 通过函数termplot() 图形化展示函数aov的分析结果. 可 以 看 出, 在R 公 式 里 面 采 用 的 运 算 符 和 表 达 式 里 面 使 用 的 运 算 符 有 着 不同的含义.例如, 公式y~x1+x2 表示模型y = β1 x1 + β2 x2 + α, 而不是(如果+ 采用它常规的含义) y = β (x1 + x2 ) + α. 为了可以在公式中使用常规的运算 符, 我们可以使用函数I: 公式y~I(x1+x2) 表示模型y = β (x1 + x2 ) + α. 类似 的, 为了定义模型y = β1 x + β2 x2 + α, 我们可以使用公式y ~ poly(x, 2) (而 非y ~ x + x^2). 但是, 为了对变量进行一定的转换, 可以在公式中包含一些 函数,如前面的杀虫剂效价分析的方差分析公式. 对于方差分析, aov() 用一个特别的语法规则定义随机效应. 例如, y ~ a + Error(b) 表示固定项a和随机项b的相加效应. 5.3 泛型函数 和 许 多 统 计 编 程 语 言 不 同 的 是, R函 数 将 输 入 对 象 的 属 性 作 为 输 入 参 数. 类是最应该关注的的一个属性.R统计函数常常返回一个类名与函数名相同的 对象(如, aov 返回类"aov"的对象, lm 返回类"lm"的对象). 我们用来解析结果 的函数对特定的类对象有特定的行为. 这些函数被称为泛型(generic)19 . 例 如, 最 常 用 的 解 析 统 计 分 析 结 果 的R函 数是summary. 它 可 以 用 来 显 示 较 为 细 致 的 结 果. 无 论 作 为 参 数 的对 象 可 能 是"lm" 类(线 性 模 型) 或 者"aov" 类(方差分析), 显示的信息显然是不一样的. 泛型函数的优势在于一个函数对 所有类的使用格式都是一样的20 . 19 译者注:在Java, C++等面向对象语言中, 泛型有更为详细的介绍.这里, 我是借用了它们的 概念.此外, 我觉得R里面的泛型,更像Java里面的接口. 20 译者注:这里和Java的接口定义非常的相似. 59 一个包含分析结果的对象常常是一个列表对象, 而它的结果展示方式由它 的类定义所决定. 前面的例子已经体现这种思想,就是一个函数的行为由输 入参数的对象类型决定21 . 这是R 22 的一个重要性质. 下面的列表列出一些用 于提取分析结果对象的信息的主要泛型函数. 这些函数的典型使用方式为: > mod <- lm(y ~ x) > df.residual(mod) [1] 8 print summary df.residual coef residuals deviance fitted logLik AIC 返回简单的汇总信息 返回较为详细的汇总信息 返回残差的自由度 返回被估计的系数(有时还包括他们的标准差) 返回残差 返回方差 返回拟合值 计算对数似然值和返回参数数目 计算Akaike 信息准则(Akaike information criterion, AIC)(依赖于logLik()) 像aov 或者lm之类的函数返回一个保存各分析结果的列表. 以前面对数据 集InsectSprays进行的方差分析为例, 我们可以看一下aov 返回对象的结构: > str(aov.spray, max.level = -1) List of 13 - attr(*, "class")= chr [1:2] "aov" "lm" 另外一种查看对象结构的方法是显示列表对象各个元素的名字: > names(aov.spray) [1] "coefficients" [4] "rank" [7] "qr" [10] "xlevels" [13] "model" "residuals" "fitted.values" "df.residual" "call" "effects" "assign" "contrasts" "terms" 这些看到的元素都可以提取出来: > aov.spray$coefficients (Intercept) sprayB 21 22 sprayC 译者注:函数多态性. 在R 里面有超过100个泛型函数. 60 sprayD 3.7606784 sprayE -1.9512174 0.1159530 sprayF 0.2579388 -2.5158217 -1.5963245 在aov()例子中的summary()函数实际上也创建了一个简单的列表对象: > str(summary(aov.spray)) List of 1 $ :Classes anova and ‘data.frame’: 2 obs. of 5 variables: ..$ Df : num [1:2] 5 66 ..$ Sum Sq : num [1:2] 88.4 26.1 ..$ Mean Sq: num [1:2] 17.688 0.395 ..$ F value: num [1:2] 44.8 NA ..$ Pr(>F) : num [1:2] 0 NA - attr(*, "class")= chr [1:2] "summary.aov" "listof" > names(summary(aov.spray)) NULL 泛 型 函 数 很 少 对对 象 进 行 任 何 操 作: 它 们 调 用 自 变 量 所 属 类 的对 应 函 数. 在R的术语里面, 泛型调用的函数称为方法(method). 简略地说, 一个方法的构 建方式是generic.cls , 其中cls 是对象的类.例如,以summary 为例, 我们可以 显示对应的方法: > apropos("^summary") [1] "summary" [3] "summary.aovlist" [5] "summary.data.frame" [7] "summary.factor" [9] "summary.glm.null" [11] "summary.lm" [13] "summary.manova" [15] "summary.mlm" [17] "summary.POSIXct" [19] "summary.table" "summary.aov" "summary.connection" "summary.default" "summary.glm" "summary.infl" "summary.lm.null" "summary.matrix" "summary.packageStatus" "summary.POSIXlt" 通过下面的简单的例子, 我们可以发现这些泛型函数在线性回归和方差分 析中的行为是不同的: > x <- y <- rnorm(5); > lm.spray <- lm(y ~ x) > names(lm.spray) [1] "coefficients" "residuals" 61 "effects" [4] "rank" "fitted.values" [7] "qr" "df.residual" [10] "call" "terms" > names(summary(lm.spray)) [1] "call" "terms" [4] "coefficients" "sigma" [7] "r.squared" "adj.r.squared" [10] "cov.unscaled" "assign" "xlevels" "model" "residuals" "df" "fstatistic" 下面的表格中显示了一些可以对分析结果对象做一些补充分析的泛型函 数, 主要参数一般都是分析结果对象, 但是有些情况下,如泛型函数如predict 或update 需要一些额外的参数. add1 drop1 step anova predict update 连续测试所有可以加入模型的元素项 连续测试所有可以从模型中移除的元素项 通过AIC (调用add1 和drop1)选择一个模型 计算一个或多个模型的方差/残差分析表 通过拟合的模型计算一个新的数据集的预测值 用新的数据或者公式拟合一个模型 还有很多各种各样用于从模型对象或者公式中提取信息的效用函数, 如函 数alias 可以用来查找一个特定公式拟合的线性模型中的线性依赖项. 最后,还有许多图形函数, 如plot 显示各种各样的诊断图, 或者termplot (见上面的例子), 尽管后面一个函数不是泛型函数但它调用泛型函数predict. 5.4 包 下表将列出R基本安装环境发布的标准 包. 其中一些包在R启动时就直 接调入内存; 这些可用通过函数search 显示: > search() [1] ".GlobalEnv" [3] "package:stats" [5] "package:grDevices" [7] "package:datasets" [9] "package:base" "package:methods" "package:graphics" "package:utils" "Autoloads" 其他包需要在载入后才能使用: > library(grid) 一个包中可以使用的函数可以通过下面的方式显示: 62 > library(help = grid) 或者直接访问网页格式的帮助文档. 任何函数的相关信息可以用( 7)页描述的 方法访问. 包 描述 base datasets grDevices graphics grid methods splines stats stats4 tcltk tools utils 基本R函数 基本R数据集 基本的或grid图形的设备函数 基本图形函数 grid图形 用于R对象和编程工具的方法和类的定义 样条回归函数和类 统计函数 基于S4标准定义的统计函数 R 和Tcl/Tk 图形接口元素的交互函数 包开发和管理的工具 R 工具函数 许 多 捐献 包 都 丰 富 了R的 统 计 方法 列 表. 它 们 各 自 发 布 , 需 要 安 装 和 载 入R. 可以在CRAN 网站上23 得到捐献包的完整列表以及相关描述. 其中的一 些包标明是推荐使用的, 因为它们包括了一些在数据分析中常用的统计方法. 推荐包常常和R的基本安装环境捆绑发布. 在下面的表中将会对它们进行简单 的描述. 23 http://cran.r-project.org/src/contrib/PACKAGES.html 63 包 描述 boot class cluster foreign 抽样和bootstraping 方法 分类方法 聚类方法 读 取 各 种 格 式(S3, Stata, SAS, Minitab, SPSS, Epi Info)的 外部数据 核密度拟合方法(包括双变量核) grid图 Venables & Ripley 著的“Modern Applied Statistics with S” 中的配套库, 包含很多有用的函数,工具和数据集 广义的可加模型 线性和非线性混合效应模型 神经网络和多项对数线性模型 递归分割 空间分析(“kriging”, 空间协方差, . . . ) 生存分析 KernSmooth lattice MASS mgcv nlme nnet rpart spatial survival 还 有 另 外 两 个 重 要 的R资 源 库: Omegahat项 目 着 力 于 网 络 应用 程 序 的 统 计计算24 并且定义R语言和一般编程语言之间的接口, Bioconductor 项目25 专 注于生物信息类的软件(特别是生物芯片数据的分析). 一个包的安装方式决定于操作系统以及R是直接从源码编译安装还是用编 译好的二进制代码安装. 在后面一种情况下,推荐采用安装CRAN 网站已经编 译好的包. 在Windows 系统中, 程序Rgui.exe 有一个“Packages”菜单, 它允许 通过网络或者本地的压缩文件直接安装. 如 果R是 本 地 编 译 的, 包 可 以 通 过 源 码(常常 是 以‘.tar.gz’作 为 扩 展 名 的 文 件) 编 译 安 装. 例 如 , 如 果 我 们 想 安 装 包gee, 首 先下 载 文 件gee 4.13-6.tar.gz (数 字4.13-6 表 明 包 的 版 本 号; 在CRAN, 一 个 包 常常 只 有 一 个 版 本 号). 然 后, 在系统控制台(不是R控制台) 键入如下命令: R CMD INSTALL gee_4.13-6.tar.gz 有几个非常有用的包管理函数, 如installed.packages, CRAN.packages, 或者download.packages. 定期键入如下命令是非常有用的: > update.packages() 这个命令将会比较系统安装包的版本与目前CRAN可以获得的包的版本(在 Windows 系统中, 这个命令可以从“Packages” 菜单中调用). 用户可以用这个 命令更新自己电脑上已经安装的包. 24 25 http://www.omegahat.org/R/ http://www.bioconductor.org/ 64 6 R编 程 实 践 自此, 我们已经对R 的功能有个全面的了解, 下面将从统计语言和编程角 度来说明. 我们将会了解一些在R编程实践中采用的简单的概念. 6.1 循环和向量化 相比下拉菜单式的程序, R的一个优势在于它可以把一系列连续的操作简 单的程序化. 这一点上和所有其他计算机编程语言是一致的, 但R有一些特性 使得非专业人士也可以很简单的编写程序. 和其他编程语言一样, R 有一些C语言类似的控制结构. 假定我们有一 个向量x, 对于向量x 中值为b的元素, 我们我们把0赋给另外一个等长向量y的 对应元素, 否则赋1. 我们首先创建一个与向量x等长的向量y: y <- numeric(length(x)) for (i in 1:length(x)) if (x[i] == b) y[i] <- 0 else y[i] <- 1 几个指令可以放在一个大括弧里面: for (i in 1:length(x)) { y[i] <- 0 ... } if (x[i] == b) { y[i] <- 0 ... } 另外一种可能的情况是当条件为真的时候一条指令才执行: while (myfun > minimum) { ... } 但 是, 由于R的 一 些 特 性, 在 很 多 情 况 下 循 环 和 控 制 结 构 可 以 避 免: 向量 化. 向量化使得循环暗含在表达式中, 前面的例子中已经多次采用了. 比如, 两 个向量之和: 65 > z <- x + y 这种加和可以通过循环结构来实现, 就像很多编程语言采用的策略一样: > z <- numeric(length(x)) > for (i in 1:length(z)) z[i] <- x[i] + y[i] 在这个例子中, 由于要用到向量的下标系统, 有必要预先创建一个向量z. 显然, 这种显式的循环仅仅用于向量x 和y 等长的情况: 如果不是这样程序代 码需要改变, 而第一种表达方式在任何情况下都成立26 . 条件语句(if ... else) 可以用逻辑索引向量代替; 同样上面的例子: > y[x == b] <- 0 > y[x != b] <- 1 除了让代码更简洁, 向量化的表达式在计算上效率更高, 特别是大数据量 的数据集. 有多个‘apply’ 形式的函数用于避免使用代码的显式循环结构. apply 作 用于矩阵的行或者/和列, 它的语法规则是apply(X, MARGIN, FUN, ...), 其 中X是一个矩阵, MARGIN 表明是对行(1) 还是列(2), 或者二者(c(1, 2))进行操 作, FUN 是 一 个 函 数(或 操 作 符, 但 是 这 种 情 况 下 操 作 符 要 在 一 个 括 弧 里 面 指 定27 ) , ... 是函数FUN可能的参数. 一个简单的例子如下. > > > > x <- rnorm(10, -5, 0.1) y <- rnorm(10, 5, 2) X <- cbind(x, y) # 矩阵X的列名保持 "x" 和 "y" apply(X, 2, mean) x y -4.975132 4.932979 > apply(X, 2, sd) x y 0.0755153 2.1388071 lapply() 可以用于一个列表对象: 它的语法类似apply 并且返回一个列 表对象. > forms <- list(y ~ x, y ~ poly(x, 2)) > lapply(forms, lm) [[1]] 26 译者注:如果对R的向量循环使用方式不了解的话,这里要小心一点. 译者注:我没有用过这种情况,原文为“or an operator, but in this case it must be specified within brackets” 27 66 Call: FUN(formula = X[[1]]) Coefficients: (Intercept) 31.683 x 5.377 [[2]] Call: FUN(formula = X[[2]]) Coefficients: (Intercept) poly(x, 2)1 4.9330 1.2181 poly(x, 2)2 -0.6037 sapply() 是函数lapply() 一个更为灵活的变种, 它可以接受向量或者矩 阵作为主要参数, 并且返回形式更为友好的结果, 常常是表格方式. 6.2 用 R写 程序 一 般 情 况 下, 一 个R程 序 以ASCII 格 式 保 存,扩 展 名 为‘.R’. 如 果 一 个工 作 要重复好多次, 用R程序是一个不错的选择. 在我们的第一个例子中, 我们想对 三个不同种属的鸟绘制一样的图, 而且数据在三个不同的文件中. 我们将一步 一步的演示看R用不同的方式去完成这个简单的问题. 首先, 我们凭直觉连续键入一系列命令, 而且预先分割图形设备. layout(matrix(1:3, 3, 1)) data <- read.table("Swal.dat") plot(data$V1, data$V2, type="l") title("swallow") data <- read.table("Wren.dat") plot(data$V1, data$V2, type="l") title("wren") data <- read.table("Dunn.dat") plot(data$V1, data$V2, type="l") title("dunnock") #分割图形界面 #读入数据 #增加标题 字符‘#’ 用于在程序中添加注释行: R 会自动跳过注释行. 67 第一个程序的问题是在我们加入另外一个物种数据时, 它过长. 此外, 一 些命令多次执行, 因此它们可以放在一起,在执行的时候仅仅修改一些参数. 这 里的策略是把参数放到一个字符型的向量中去, 然后用下标去访问这些不同的 值. 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)) { data <- read.table(file[i]) # 读入数据 plot(data$V1, data$V2, type="l") title(species[i]) # 增加标题 } 注 意 代 码read.table()里 面 的 参 数file[i] 上 面没 有 双 引 号, 因 为 这 个 参数是字符型. 现在的代码比较紧凑. 它比较容易加入新的物种, 因为设置物种名字和数 据文件的向量都程序的前端. 如 果 扩 展 名 为‘.dat’的 数 据 文 件 在R的 工 作 目 录 下 面, 程 序 可 以 正 常 运 行, 否则用户要设置工作目录, 或者设置绝对路径(例如: file <- "/home/paradis/ data/Swal.dat"). 如 果 程 序 保 存 在 文 件Mybirds.R 中, 可 以 通 过 键 入如 下 命 令执行: > source("Mybirds.R") 和所有以文件作为输入对象的函数一样, 如果该文件不在当前工作目录下 面, 用户需要提供该文件的绝对路径. 6.3 编写你自己的函数 大多数R的工作是通过函数来实现, 而且这些函数的输入参数都放在一个 括弧里面. 用户可以编写他们自己的函数, 并且这些函数和R里面的其他函数 有一样的特性. 编写自己的函数可以让你有效的,灵活的与合理的使用R. 我们再次使用前 面读数据并且画图的例子. 如果我们想在其他情况下进行这样的操作, 写一个 函数是一个不错的想法: myfun <- function(S, F) { data <- read.table(F) plot(data$V1, data$V2, type="l") title(S) } 68 执 行 时, 这 个 函 数 必 须 载 入 内 存.当 然, 这 有 多 种 方 式实 现. 和 所 有 其 他 命令一样, 函数的各行可以直接通过键盘键入,或者从一个文本编辑器里面拷 贝 粘 贴.如 果函 数 保 存 在 一 个 文 本 文 件 中, 可 以 命 令source() 载 入. 如 果 用 户 期 望 一 些 函 数 在R启 动 时 就 被 载 入, 可 以 把 它 们 保 存 在 工 作 目 录 下 面 的 文 件.RData中. 另外一种方式是, 配置文件‘.Rprofile’ 或‘Rprofile’ (详见?Startup for details). 最后, 还可以创建一个包, 但是这里不想多讨论(见手册“编写R扩 展”). 一旦函数载入后, 我们就可以键入一条命令以读入数据和画出我们想要的 图, 如myfun("swallow", "Swal.dat"). 因此, 现在我们的程序有第三个实现 的版本了: layout(matrix(1:3, 3, 1)) myfun("swallow", "Swal.dat") myfun("wren", "Wrenn.dat") myfun("dunnock", "Dunn.dat") 我们还可以用sapply() 实现程序的第四个版本: layout(matrix(1:3, 3, 1)) species <- c("swallow", "wren", "dunnock") file <- c("Swal.dat" , "Wren.dat", "Dunn.dat") sapply(species, myfun, file) 在R里 面, 没 有 必 要 在 一 个 函 数 里 面 进 行 变 量 声 明. 当 一 个 函 数 执 行 时, R用 一 种 称 为 词汇作用域(lexical scoping)的 规 则 决 定 一 个 对 象 的 作 用域 相 对 一个函数是局部还是全局.为了理解这种机制, 我们可以认真研究一下下面的 一个简单函数: > foo <- function() print(x) > x <- 1 > foo() [1] 1 x 不是为了在函数foo()里面创建对象, 因此R 将会在封装环境中搜索是 否有个名为x 的对象, 和打印它的值(否则一条错误信息将会显示, 执行中断). 如果x是我们定义的函数中一个对象的名字, 全局环境中变量x 值将会被 采用. > x <- 1 > foo2 <- function() { x <- 2; print(x) } > foo2() [1] 2 >x [1] 1 69 此时, print() 使用在它所在的环境中定义的x, 即环境foo2 中的x. 前面提及的“封装”是关键所在.在前面两个演示函数中, 有两个环境: 全局 环境和函数foo 或foo2的局部环境. 如果有三个或者更多的嵌套环境, 对象搜 索将逐层搜索直到全局环境. 有两种方式指定一个函数的参数: 通过它们的定义时的位置或者名字(又 称为标签参数). 例如, 假定一个函数有三个参数: foo <- function(arg1, arg2, arg3) {...} foo() 在执行时可以不用名字arg1, . . . , 如果相应的参数对象放在相应的 位置上, 如: foo(x, y, z). 但是, 如果使用了参数的名字, 位置信息将会失效, 如foo(arg3 = z, arg2 = y, arg1 = x). R函数的另外一个特性是函数可能 采用定义时的默认设置. 例如: foo <- function(arg1, arg2 = 5, arg3 = FALSE) {...} 命 令foo(x), foo(x, 5, FALSE) 和foo(x, arg3 = FALSE) 将 会 得到 一样 的 结 果. 使 用 一 个 函 数 的 默 认 设 置 非 常 有 用, 特 别 在 使 用 标 签 参 数 的 时 候, 如foo(x, arg3 = TRUE) 仅仅改变一个默认设置. 在结束本章前, 我们来看另外一个例子. 尽管这个例子不是纯的统计学例 子, 但是它很好地展示了R 语言的灵活性. 假定我们想研究一个非线性模型的 行为: 这个模型(Ricker 模型)的定义如下: Nt+1 = Nt exp r 1 − Nt K 这 个 模 型 广 泛 地 用于 种 群 动 态 变 化 (population dynamics) 的 研 究 里 面,特 别 是 鱼 类 的 种 群 变 化. 我 们 想 用 一 个 函 数 去 模拟 这 个 模 型 关 于 增 长 率r 和初始群体大小N0 的变化情况(承载能力K 常常设定为1且以这个值作为默认 值)的影响; 结果将以种群大小相对时间的图表示. 我们还将设定一个可选项允 许用户只显示最后若干步中种群大小(默认所有结果都会被绘制出来). 下面的 函数就是做Ricker模型的数值模拟的. ricker <- function(nzero, r, K=1, time=100, from=0, to=time) { N <- numeric(time+1) N[1] <- nzero for (i in 1:time) N[i+1] <- N[i]*exp(r*(1 - N[i]/K)) Time <- 0:time plot(Time, N, type="l", xlim=c(from, to)) } 你可以试一试下面的代码: 70 > > > > layout(matrix(1:3, 3, 1)) ricker(0.1, 1); title("r = 1") ricker(0.1, 2); title("r = 2") ricker(0.1, 3); title("r = 3") 71 7 R 相关的文献 手册. 在目录R HOME/doc/manual/下面有几份R安装时分发的手册: • R导论28 [R-intro.pdf], • R 安装和管理 [R-admin.pdf], • R 数据导入/导出 [R-data.pdf], • 编写R扩展 [R-exts.pdf], • R 语言定义 [R-lang.pdf]. 这些文件可能以不同的格式(pdf, html, texi, . . . )显示, 这取决于你采用 的安装方式. FAQ. 在 目 录R HOME/doc/html/ 还 有 一 个FAQ (常见问题集, Frequently Asked Questions) .R-FAQ 的版本会在CRAN网站定期更新: http://cran.r-project.org/doc/FAQ/R-FAQ.html 在线资 源. CRAN 网站拥有许多文档, 参考文献和其他有用的R 网页链接. 它 还 有有 关R或 者 统 计 方法 的 出 版 物(书 或 者 文 章)的 清 单29 以 及 一 些R用 户编写的文档和手册30 . 邮件列 表. R有 四 种 形 式 的 邮 件 讨 论列 表; 订 购, 发 布 信 息 以 及 阅 读 所 有 档 案 记录请参阅: http://www.R-project.org/mail.html. 一 般 的 讨 论列 表 是‘r-help’.对 于R用 户 来 说, 这 是 一 个 非 常 有 意 思 的 资 源(其他三个列表主要面向开发用户和新版本发布).许多用户向‘r-help’咨 询 的 函 数 或 程 序 在 档 案 文 件 中 都 已 有 回 答. 如 果 在 使 用R的 过 程 中 碰 到 困难, 你首先应该根据下面的建议处理一下, 然后才考虑是否给‘r-help’发 咨询信息: 1. 仔细阅读在线帮助文档(可能需要搜索引擎); 2. 阅读R-FAQ; 3. 搜索‘r-help’档案文件, 或者使用一些网站提供的搜索引擎31 ; 4. 在你提交问题前阅读“公告指南”32 . 28 译者注:已经有中文版发布.用google 搜索即可得到. http://www.R-project.org/doc/bib/R-publications.html 30 http://cran.r-project.org/other-docs.html 31 这些网站的地址清单在http://cran.r-project.org/search.html 32 http://www.r-project.org/posting-guide.html 29 72 R 新 闻. 电 子 杂 志R News 的 目 的 是 弥 补 电 子 讨 论列 表 和 传 统 科 学 文 献 的 鸿 沟. 第一期是在2001年一月发布的33 . 在文献 中引用 R. 最后, 如果你想在你的文章中引用R, 你必须采用下面参考文 献: R Development Core Team (2005). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL: http://www.R-project.org. 33 http://cran.r-project.org/doc/Rnews/ 73 ...
View Full Document

Ask a homework question - tutors are online