R语言基础
写在前边
自2019年9月30日以来，本页面已被访问20004人次。
为了方便学习，教程中把执行代码后打印在屏幕上的信息，附在代码后边，前边加上两个#号及两个空格。练习时不用输入#号及其后边的内容。比如：
print("欢迎一起学习R语言！") ## [1] "欢迎一起学习R语言！"
这里[1] “欢迎一起学习R语言！”就是print(“欢迎一起学习R语言！”)这条命令执行后打印到屏幕上的信息。
R语言简介
若对R语言的历史和特点等感兴趣可以看看R官网的介绍。
R的安装
这里只讲Windows下的安装，能用Linux系统的基本上也就不太可能不会安装软件。我没用过Mac，不明确其中可能会遇到的问题，这里就不介绍了。
到R官网或者国内（比如清华或同济）的镜像下载安装文件进行安装。建议不要安在C盘，不然可能遇到一些权限问题，安装路径中避免有非英文字符。如果已经安装在C盘，当安装R包的时候如果遇到不可写的问题，则需要修改安装目录下library文件夹的权限。 启动选项选择“Yes(自定义启动)”，显示模式选择“SDI（多个分开的窗户）”，帮助风格选择“HTML帮助”，这样的选择以后用起来可能会方便一些。可以观看R的安装微课进一步熟悉。
R包的安装详见R包安装一节。
基本装备
通常，为了在使用R的过程中更轻松和愉快，可以添加一些装备，我喜欢用独立的纯文本编辑器。Notepad++作为纯文本编辑器与NppToR组合，对于编写代码和推送到R控制台执行都是非常方便的。Notepad++在GitHub上的页面有时候连不上，可以试试Notepad++在SourceForge上的镜像和NppToR在SourceForge上的镜像。可以观看R的安装微课进一步熟悉。
R控制台
控制台的上边的菜单栏里有很多的选项，自己去尝试点开看看都有哪些功能。比如尝试一下把控制台的字号改大一点。(注：这里是个示例，并不是这个4.2.0版有什么特殊的。大家安装最新的版本用于学习即可。)
如上图所示，在R控制台中，>后边输入新的命令，当>出现时，说明先前的命令已经执行完成，等待新命令的输入。命令输入完成之后，按Enter键命令就得到了执行，当命令执行完成后，出现新的>。
而+出现时，说明当前的命令输入不完整，这时候如果按Enter键将会不断出现+，直到命令输入完整后，按Enter键才执行。 当自认为命令输入完整后，依然出现+时，可能是什么地方漏掉了括号、引号等成对出现的符号，从而导致控制台认为输入不完整，应该细心检查。
命令执行完成后，有的会在屏幕上打印一些信息，这些信息根据不同的命令而有所不同。对于打印出的结果是向量的时候，每一行会显示本行第一个元素的索引。比如上图中1:100这个命令创建了一个包含1~100这100个数字的向量，这个命令运算结果并没有赋值给任何对象，默认就是将它print到屏幕上，因此我们看到print出来的结果。
当命令中出现错误的时候，错误信息会展示出来，另外还可能有警告信息。
R控制台里有一些快捷键可以使用：
| Esc | 退出当前正在执行的进程，当陷入死循环或者代码输入有错一直显示未输入完成时可用此键退出 |
| Ctrl + X | 将选择的代码复制并粘贴到输入位置 |
| Ctrl + L | 清屏， 不再显示以往输入的代码和打印的结果，不影响先前定义的对象 |
| ↑ | 翻到先前执行过的命令 |
| ↓ | 向下翻后续执行过的命令 |
| Ctrl + A | 光标转到行首 |
| Ctrl + E | 光标转到行尾 |
寻求帮助
R内
通过help序列的函数，可以调出R内在的对象的文档。在安装的时候选取“HTML帮助”的话，这个帮助手册还有许多超链接到相关的内容。
? c ? t.test
?后边接函数名，用于查看其使用说明，比如这里? c 查看c函数的使用说明。? t.test查看t.test函数的使用手册，这等同于help(t.test)，
? "for" ? "+"
对于操作符以及for、if、function等具有句法意义的词，必须要加上引号，不然可能报错或者后边出现+号说明命令未输入完成。如果报错，就试试加上引号。
?? plot
当查询的函数所在的包没有加载，可以用??来搜索，这等同于help.search(“plot”), help.search里边的引号是必须加上的。
有的函数的用法可以用example函数调出示例，比如example(plot)可以查看plot函数的示例，这里注意留意R控制台的提示和画图设备上的效果。
help.start()命令，可以调出一个页面，此页面上有连接到别的帮助页面。 其中“An Introduction to R”是很好的入门指南，而“Search Engine & Keywords”也很有用。 开卷有益，页面上的连接都可以点开看看。
RSiteSearch函数用于在R网站上搜索，比如RSiteSearch(“ggplot”)可搜索ggplot这个函数，注意，这里的引号不能省略。
在往后的学习或工作中遇到不理解的地方，即可使用这节的方法快速查看参考文档。
R外
用好搜索引擎，很多时候我们所遇到的问题都是别人遇到过的，网上搜一下很容易就找到答案了。 也可以到StackOverflow上进行搜索或者提问。
另外一个重要的提问的地方是我们的QQ，每个学期会有一个QQ群，大家可以相互交流。
命令构建
R的基本命令由表达式、赋值之一或者两者组成：单独的表达式作为命令的时候，表达式将被执行并且将结果打印到屏幕上（除非专门设置不打印），也就是调用了print函数，所获得的结果将丢失；而赋值命令将执行表达式，并且将所得的值传递给一个对象，但是结果不会自动打印（除非专门设置）。
R命令中的字母是区分大小写的，字符是区分全角和半角的，输入的时候要小心，避免输错。当命令输错的时候，会有报错信息，仔细阅读，看看错误的地方在哪里。
3 + 5 - 2 * 0 ## [1] 8
这是一个表达式，作为命令时将被执行，然后结果print到屏幕上，然后这个计算结果没法调用了。
x <- 3 + 5 - 2 * 0 x ## [1] 8 x^2 ## [1] 64 (y <- sin(pi*3/4)) ## [1] 0.7071068 y ## [1] 0.7071068
表达式运算的结果被赋值给一个叫x的对象，不再将其值print到屏幕上，但是结果可以通过x用于后续的运算 。 如果要将运算结果的信息print的屏幕上，可以将x作为命令。 如果要同时复制且又把信息print到屏幕上，那么可以把命令用小括号括起来。
对于不同class的对象，print的信息是不一样的，有的print整个对象的值，比如这里对于一个向量，就是把向量所有值print出来；有的print的是一些特定的信息，这个在后边的学习中会进一步体会。
多个命令可以放在同一行，中间用分号;隔开，将依次执行：
x <- 1:10; y <- 2:5; x; y^2 ## [1] 1 2 3 4 5 6 7 8 9 10 ## [1] 4 9 16 25
多行命令要在同等条件下执行，使用花括号{ } 括起来：
{ c(1,2,3) 9:7 5:7 } ## [1] 5 6 7
花括号用于将一组命令组合在一起，可以看到只有最后一个命名的结果会被值被print出来了。
{ print(c(1,2,3)) 9:7 5:7 } ## [1] 1 2 3 ## [1] 5 6 7
print不管放在哪里都会把结果print到屏幕。
R里赋值习惯使用←，而使用=和→也可以的，另外还能连续赋值。试试：
Y = 1:5 Y ## [1] 1 2 3 4 5 1:5 -> Z Z ## [1] 1 2 3 4 5 X <- Y <- Z <- 1:10 X;Y;Z ## [1] 1 2 3 4 5 6 7 8 9 10 ## [1] 1 2 3 4 5 6 7 8 9 10 ## [1] 1 2 3 4 5 6 7 8 9 10
除了用赋值运算符以外，还可以用assign函数，如下：
assign("object_name", expression)
这里的“object_name”是用引号引起来的,这里需要一个字符，这和不引的时候是有区别的，这就牵涉到R的数据类型（data type）的问题。
R的对象
所有的计算机语言都通过变量（variable）取用内存中的数据，R中的变量不直接取用内存中的数据，而是取用专门数据结构——称为对象（object）。R里的一切都是对象，对象通过符号（symbols）或变量指及，R中发生的一切都是函数调用（function call）。函数(function)本身也是对象。 function也可以是运算符的形式，也可以省略（也就是直接在控制台输入对象或者其名称时，就会返回相应的结果，比如对象的值或者摘要）。这样在编程的时候实际上是处理对象之间的相互作用，称为面向对象的编程（Object Oriented Programming, OOP）。
对象的模式（mode）和长度（length）
每个对象都有模式（mode）和长度（length）这两个内在属性(Intrinsic attributes)。对象x模式用mode(x)查看，长度用length(x)查看。
模式(Mode)是基于对象的类型（Type）的。对象x的类型用typeof(x)函数查看，是对象在R语言内部的类型及储存模式。而mode()则是根据typeof()进行的互斥的归类。在使用者层面通常使用mode(), 而不去关注R语言层面更细分的typeof()。每个对象只属于一个mode。在typeof层面区分整数（integer）、实数（double），而mode层面都是数值型(numeric)。mode层面的函数（function）在typeof层面又区分为closure、builtin、special等类型。很显然，我们在后续在讨论对象的类型的时候，是在mode层面来说的。
x <- 1:10 x ## [1] 1 2 3 4 5 6 7 8 9 10 mode(x) ## [1] "numeric" length(x) ## [1] 10 y <- character(length=0) y ## character(0) length(y) ## [1] 0 mode(y) ## [1] "character" length(x) <- 20 x ## [1] 1 2 3 4 5 6 7 8 9 10 NA NA NA NA NA NA NA NA NA NA length(x) <- 5 x ## [1] 1 2 3 4 5 mode(x) <- "character" x ## [1] "1" "2" "3" "4" "5" typeof(x) ## [1] "character" typeof(x) <- "integer" ## Error in typeof(x) <- "integer" : 没有"typeof<-"这个函数 mode(x) <- "integer" mode(x) ## [1] "numeric"
对象的类别 (class)
每个对象都至少属于一个类别（class），class可以有多个，针对每种class定义了一些方法，也就是函数的执行形式。对于一些一般的函数，用于不同的类的对象效果不同，这是在定义类的时候就定义了这个函数的用法。 具有同样class的对象在一些操作上可以同等对待。 对于简单的向量，类别(class)就是其模式(mode)，元素全为整数的向量其class为integer，而mode依然是numeric。
x <- 1:10 class(x) ## [1] "integer" set.seed(2022); y <- sample(1:10) y ## [1] 4 3 6 7 9 10 2 8 5 1 fit <- lm(y~x) class(fit) ## [1] "lm"
print(x) ## [1] 1 2 3 4 5 6 7 8 9 10 print(fit) ## ## Call: ## lm(formula = y ~ x) ## ## Coefficients: ## (Intercept) x ## 4.7333 0.1394 ##
print函数对不同的class的对象效果不一样的。
x <- 1:10 class(x) ## [1] "integer" class(x) <- "lm" class(x) ## [1] "lm" str(x) ## 'lm' int [1:10] 1 2 3 4 5 6 7 8 9 10 print(x) ## 错误: $ operator is invalid for atomic vectors x ## 错误: $ operator is invalid for atomic vectors
改变对象的class并没有改变对象的结构，print按照class标签进行而数据结构不对，就不能够正确print。
X <- data.frame(A=1:10, B=y) mode(X) ## [1] "list" class(X) ## [1] "data.frame"
数据框的mode是列表，而class是数据框。
对象的附加属性(attributes)
一些对象还具有非内在的附加属性(attributes)，也就是给对象加上的一些标签，class也是一种attribute。通过attributes函数可以查看或添加，通过attr函数可以增加特定的attributes。通过attributes和attr查看和改变的只是属性的标签，这个标签申明的属性可以不是对象本身所具备或者适用的。正常情况下对象具有相应的属性才会有相应的标签。
x <- 1:10 attributes(x) ## NULL attr(x, "dim") ## NULL
对于简单的向量，没有其它的属性，class用class()函数查看。
attr(x, "dim") <- c(2, 5) x ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 class(x) ## [1] "matrix" "array" attributes(x) <- NULL x ## [1] 1 2 3 4 5 6 7 8 9 10 class(x) ## [1] "integer"
添加了dim这个属性，dim即dimension，这里指定x为2行5列。增加dim属性之后，x的class也变成了matrix，这时候，适合于matrix的方法将适合于x 。而去掉了sttributes，x又变回了简单的向量。
X <- data.frame(A = 1:10, B=letters[1:10]) mode(X) ## [1] "list" class(X) ## [1] "data.frame" attributes(X) ## $names ## [1] "A" "B" ## ## $class ## [1] "data.frame" ## ## $row.names ## [1] 1 2 3 4 5 6 7 8 9 10 attributes(X) <- NULL X ## [[1]] ## [1] 1 2 3 4 5 6 7 8 9 10 ## ## [[2]] ## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" class(X) ## [1] "list"
可以看到，数据框的mode是list，而其attributes包括names、class、row.names这些，如果去掉了这些attributes，就变成了一般的list。
(x <- 1:10) attributes(x) <- list(names="A", row.names=c("B", "C"), class="data.frame") mode(x) ## [1] "numeric" class(x) ## [1] "data.frame" print(x) ## > x ## A NA NA NA NA NA NA NA NA NA ## B 1 2 3 4 5 6 7 8 9 10 ## C <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> ## Warning message: ## In format.data.frame(if (omit) x[seq_len(n0), , drop = FALSE] else x, : ## 破损的数据框：縮短列的长度或用NA来填充空白
而如果强行的给一个向量data.frame的这些attributes，那么并不会改变它的mode，而其class发生了改变，这也改变了print的效果。
对象的结构
对象都有内部结构structure，通过str函数查看
x <- 1:10 str(x) ## int [1:10] 1 2 3 4 5 6 7 8 9 10 X <- data.frame(A = 1:10, B=letters[1:10]) str(X) ## 'data.frame': 10 obs. of 2 variables: ## $ A: int 1 2 3 4 5 6 7 8 9 10 ## $ B: chr "a" "b" "c" "d" ...
需要注意，structure函数是用来根据附加的属性来结构化数据的，并不是用来查看数据的结构的。
x <- 1:10 y <- structure(x, dim=c(2, 5)) y ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 str(y) ## int [1:2, 1:5] 1 2 3 4 5 6 7 8 9 10 z <- str(x) ## int [1:10] 1 2 3 4 5 6 7 8 9 10 str(z) ## NULL
str函数把结果打印到屏幕上，而不会返回一个值，所以我们看到并没有值赋值给z。
值类型
mode
R的值类型（value type）有数值(Numeric)、字符(Character) 、复数（Complex）逻辑（Logical） 等。R里并没有标量（Scalar），而单一的值被看作长度为1的向量。向量是原子型的，向量的类型就是其中值的类型。使用mode查看向量中的值类型。不同类型的值能够进行的运算或者操作是不一样的，比如字符型的就没法求和，也没法开平方。
(x <- c(1, 3, 5)) ## [1] 1 3 5 mode(x) ## [1] "numeric" is.numeric(x) ## [1] TRUE sum(x) ## [1] 9
用is.numeric来判断x是否为数值型，同样可以使用is.character、is.complex、is.logical等判断值是否是字符、复数、逻辑型。
x是数值型的，可以进行针对数值型的运算，这里x所有元素的加和,类似的mean、median、var、sd函数秋x的均值、中位数、方差、标准差。
(y <- as.character(x)) ## [1] "1" "3" "5" sum(y) ## Error in sum(y) : 'type'(character)参数不对
字符型的值都是带引号的。使用as.character将向量x变为数值型，原先x并没有改变，而是新产生一个向量，如果不赋值给x，就不会给变x的值。同样也可以使用as.numeric、as.complex、as.logical等等将其它类型的值转换为数值、复数、逻辑型的值。
字符型的向量不能使用针对数值型的函数，因为没有办法对字符型的向量求和（比如）。
(X <- c("1", "0.001", "A", NA)) ## [1] "1" "0.001" "A" NA as.numeric(X) ## [1] 1.000 0.001 NA NA ## Warning message: ## 强制改变过程中产生了NA as.complex(X) ## [1] 1.000+0i 0.001+0i NA NA ## Warning message: ## 强制改变过程中产生了NA
字符型的向量转为数值型的时候，如果字符不是数字的形式，那么没有对应的数字，转化的时候就会产生NA。
(z <- c(-1, 0, 1, 0.0003, NA)) ## [1] -1.0000 0.0000 1.0000 0.0003 NA (Z <- as.logical(z)) ## [1] TRUE FALSE TRUE TRUE NA as.numeric(Z) ## [1] 1 0 1 1 NA as.character(Z) ## [1] "TRUE" "FALSE" "TRUE" "TRUE" NA
在数值转化为逻辑型的时候，0是FALSE，非0都是TRUE。 逻辑型只有TRUE和FALSE两个值。逻辑值也存在NA值。转化为数值时，变成0和1，NA依然是NA。 转换成字符型则加上了引号，NA除外。
sqrt(-1) ## [1] NaN ## Warning message: ## In sqrt(-1) : 产生了NaNs sqrt(as.complex(-1)) ## [1] 0+1i sqrt(-1 + 0i) ## [1] 0+1i 1 + 2 ## [1] 3 1 + "2" ## Error in 1 + "2" : 二进列运算符中有非数值参数 TRUE + 0 ## [1] 1 TRUE & 1i ## [1] TRUE TRUE & "1" ## Error in TRUE & "1" : 此运算只适用于数值，逻辑值或复数值种类
数值型的-1是没法开根号的，转成复数型的，就可以开根号了。+ xi （x为数值）这样的形式，也能得到复数型的值。 数值型之间可以相加等数学运算，字符型的不能进行数学运算。逻辑型的可以和数值型进行数学运算，逻辑型的可以和数值型、复数型的进行逻辑运算，字符型的不能进行逻辑运算。这里不一一列举，如果有不能进行的运算，要留意报错信息。
(x <- c(TRUE, FALSE, 0)) ## [1] 1 0 0 mode(x) ## [1] "numeric" (y <- c(TRUE, FALSE, 0, 10 + 0i)) ## [1] 1+0i 0+0i 0+0i 10+0i mode(y) ## [1] "complex" (z <- c(TRUE, FALSE, 0, 10 + 0i, "A")) ## [1] "TRUE" "FALSE" "0" "10+0i" "A" mode(z) ## [1] "character"
在使用c()函数创建向量的时候，如果将不同类型的值放在一起，会有一个转化，这种优先级是这样的：字符型 > 复数型 > 数值型 > 逻辑型 。一个向量只有一种类型，在列表读入的时候，也会遇到这种同一列里可能有几种类型的情况，这时候也会有转化，比如一列数字里如果具有含有标点的字符，这时候就会被当作字符型的读入。 而在老版本的R中，默认情情况下会将字符转为因子型，这样如果直接使用as.numeric转换成数值型，就会发生错误，需要先转成字符型再转成数字型。现在的新版本R中不存在这个问题了，大家在一些老机器上如果遇到一些旧版本的R需要留意。
type
除了mode和class以外，对象的类型还有typeof和storage.mode两个函数返回的值，typeof函数和mode函数的类型大都具有相同的名称，除了：typeof的integer和double在mode里都是numeric；typeof的special和builtin在mode里是function；typeof里的symbol在mode里叫name；typeof的language在mode里是(或者call。storage.mode通常和mode的值是一样的，而是这两个都是基于typeof的。typeof是R内部的type，storage.mode是mode初级版本，更高效，也更接近于底层，比如storage.mode的integer和double和typeof是一样的。
typeof ## function (x) ## .Internal(typeof(x)) ## <bytecode: 0x00000257fac88030> ## <environment: namespace:base>
可以看到typeof函数是一个内部函数，具体怎么实现的就不是我们这样的运用层面需要去弄清楚的了。
mode ## function (x) ## { ## if (is.expression(x)) ## return("expression") ## if (is.call(x)) ## return(switch(deparse(x[[1L]])[1L], `(` = "(", "call")) ## if (is.name(x)) ## "name" ## else switch(tx <- typeof(x), double = , integer = "numeric", ## closure = , builtin = , special = "function", tx) ## } ## <bytecode: 0x00000257fabe4c38> ## <environment: namespace:base>
可以看到mode是基于typeof的。
storage.mode ## function (x) ## switch(tx <- typeof(x), closure = , builtin = , special = "function", ## tx) ## <bytecode: 0x00000257fae5b978> ## <environment: namespace:base>
storage.mode 函数也是基于typeof函数返回的结果的。
一些特殊的值（常量）
NA
NA意为Not Available，表示变量值未知。数据中包含缺失的值的时候，R将之设置为NA，比如将一个不是数字样式的字符强行转为数值型的，没有与之对应的数字，得到NA。NA不等于任何值。NA默认是一个逻辑型的。
as.numeric(c("3.14", "pi")) ## [1] 3.14 NA ## Warning message: ## 强制改变过程中产生了NA mode(NA) ## [1] "logical" 1==NA ## [1] NA TRUE & NA ## [1] NA TRUE | NA ## [1] TRUE FALSE | NA ## [1] NA is.na(NA) ## [1] TRUE 1 + NA ## [1] NA
通过is.na来判断一个值是否是NA，NA和数值型的进行运算，得到的还是NA，需要留意的是TRUE | NA得到的结果是TRUE。
NaN
NaN意为Not a Number，表示未定义的数学运算结果。
Inf - Inf ## [1] NaN 0/0 ## [1] NaN mode(NaN) ## [1] "numeric" is.na(NaN) ## [1] TRUE
NaN是NA的一种特殊形式。
NULL
NULL是R里的一个特殊对象，其值为NULL，由NULL指向，其模式为NULL， 用于表示一个对象或者元素是不存在的。
NULL ## NULL NULL <- 0 ## Error in NULL <- 0 : (do_set)赋值公式左手不对 is.null(NULL) ## [1] TRUE is.na(NULL) ## logical(0) c(1, 2, NULL) ## [1] 1 2 c(1, 2, NA) ## [1] 1 2 NA
没法给NULL赋值，NULL是一个特殊的保留的符号。
Inf
计算结果超过R数字的极限时，结果为Inf，意为Infinity。通常因为除以0而产生。
10^1000 ## [1] Inf 1/0 ## [1] Inf -1/0 ## [1] -Inf
对象类型 (I)
R的对象有许多类型，其中向量（vector）称为原子型的结构，因为同一个向量中的所有元素都只能是一个模式。向量中的元素可以是数值型(numeric)、复数型(complex)、逻辑型(logical)、字符型(character)以及原始型（raw）。
另外R中的对象类型还有因子(factor)、数据框(data frame)、矩阵(array)、列表(list)、函数(function)、表达式（expression）、调用(call)、S3、S4等等的。
向量
向量（vector）可以看作一系列装着数据的连续的小格，可以通过小格的序号来取出数据。可以通过c来创建。seq等函数的结果也是向量。
向量是原子型的结构，同一个向量中所有的元素都只能是一个类型，向量的类型就是其元素的类型。 而每个元素单独拿出来，实际上就是一个长度为1的向量。 给向量元素加上名称之后，向量可当字典用。
(x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) ## [1] 1 2 3 4 5 6 7 8 9 10 str(x) ## num [1:10] 1 2 3 4 5 6 7 8 9 10 mode(x) ## [1] "numeric" class(x) ## [1] "numeric" attributes(x) ## NULL names(x) <- LETTERS[1:10] x ## A B C D E F G H I J ## 1 2 3 4 5 6 7 8 9 10 str(x) ## Named num [1:10] 1 2 3 4 5 6 7 8 9 10 ## - attr(*, "names")= chr [1:10] "A" "B" "C" "D" ... mode(x) ## [1] "numeric" class(x) ## [1] "numeric" x["A"] ## A ## 1 x[c("C", "D", "B", "A")] ## C D B A ## 3 4 2 1 is.vector(x) ## [1] TRUE is.numeric(x) ## [1] TRUE (y <- as.character(x)) ## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" str(y) ## chr [1:10] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" mode(y) ## [1] "character" class(y) ## [1] "character" (x <- c(c("cat", "dog"), "mouse", c("human"))) ## [1] "cat" "dog" "mouse" "human"
根据数值的类型，R里有六种基本的向量类型：字符型(character)、逻辑型(logical)、整数型(integer)、实数型(double)、复数型(complex)、原始型(raw)。而整数型和实数型又都是数值型（numeric）。用is序列函数（is.character、is.numeric、is.integer、is.double、is.complex、is.raw）来判断是否是相应的类型，而用as系列函数（as.character、as.numeric、as.integer、as.double、as.complex、as.raw）来在这些类型间进行转换，转换时需要注意是否确实可以转换。
矩阵
矩阵（matrix）是按行和列排列的一组数值。 给向量加上dim()附加属性就可以得到矩阵。通过matrix函数创建。
(x <- matrix(1:10, nrow=2)) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 mode(x) ## [1] "numeric" class(x) ## [1] "matrix" "array" attributes(x) ## $dim ## [1] 2 5 str(x) ## int [1:2, 1:5] 1 2 3 4 5 6 7 8 9 10 length(x) ## [1] 10 nrow(x) ## [1] 2 ncol(x) ## [1] 5
可以使用matrix()来创建矩阵，矩阵的mode还是numeric，本质上是向量。其class为matrix和array ，具有一般向量所不具有的特征。
x <- 1:10 dim(x) <- c(2, 5) x ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 str(x) ## int [1:2, 1:5] 1 2 3 4 5 6 7 8 9 10 rownames(x) <- letters[1:2] colnames(x) <- LETTERS[1:5] x ## A B C D E ## a 1 3 5 7 9 ## b 2 4 6 8 10 attributes(x) ## $dim ## [1] 2 5 ## ## $dimnames ## $dimnames[[1]] ## [1] "a" "b" ## ## $dimnames[[2]] ## [1] "A" "B" "C" "D" "E"
通过给一个向量dim个attributes也可以创建矩阵。通过dim函数来设置dim附加属性，也可以用attr(x, “dim”) ← c(2, 5)，另外，有dimnames这一attributes，可以通过rownames和colnames来获取或者添加。 。
(x <- matrix(1:10, nrow=2, byrow=TRUE)) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 2 3 4 5 ## [2,] 6 7 8 9 10 attributes(x) <- NULL x ## [1] 1 6 2 7 3 8 4 9 5 10
去掉矩阵的attributes之后，矩阵变成了一般的向量，还是先按列排，第一列然后第二列，以此类推。使用matrix函数创建矩阵的时候，byrow=TRUE使得元素的排列先排行再排列，在创建matrix的时候byrow=TRUE只是在创建的时候起效，创建之后这一信息并未保存在对象信息中。 # 显然，
数组
矩阵可以看作二维的数组(array) ，因此推广一下理解更高维数组。在平时的工作中很少用到。
(x <- array(1:10, c(2, 5))) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 (x <- array(1:24, c(2, 3, 4))) ## , , 1 ## ## [,1] [,2] [,3] ## [1,] 1 3 5 ## [2,] 2 4 6 ## ## , , 2 ## ## [,1] [,2] [,3] ## [1,] 7 9 11 ## [2,] 8 10 12 ## ## , , 3 ## ## [,1] [,2] [,3] ## [1,] 13 15 17 ## [2,] 14 16 18 ## ## , , 4 ## ## [,1] [,2] [,3] ## [1,] 19 21 23 ## [2,] 20 22 24 (x <- array(1:40, c(2, 5, 2, 2))) ## , , 1, 1 ## ## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 3 5 7 9 ## [2,] 2 4 6 8 10 ## ## , , 2, 1 ## ## [,1] [,2] [,3] [,4] [,5] ## [1,] 11 13 15 17 19 ## [2,] 12 14 16 18 20 ## ## , , 1, 2 ## ## [,1] [,2] [,3] [,4] [,5] ## [1,] 21 23 25 27 29 ## [2,] 22 24 26 28 30 ## ## , , 2, 2 ## ## [,1] [,2] [,3] [,4] [,5] ## [1,] 31 33 35 37 39 ## [2,] 32 34 36 38 40 mode(x) ## [1] "numeric" class(x) ## [1] "array" attributes(x) ## $dim ## [1] 2 5 2 2 str(x) ## int [1:2, 1:5, 1:2, 1:2] 1 2 3 4 5 6 7 8 9 10 ...
列表
列表(list) 可以看作一般化的向量，只是其元素可以是不同类型，甚至是list。
(x <- list(A=1:3, B=letters[1:10]) ) ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" str(x) ## List of 2 ## $ A: int [1:3] 1 2 3 ## $ B: chr [1:10] "a" "b" "c" "d" ... mode(x) ## [1] "list" class(x) ## [1] "list" attributes(x) ## $names ## [1] "A" "B" attributes(x) <- NULL x ## [[1]] ## [1] 1 2 3 ## ## [[2]] ## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j"
数据框
数据框(data frame)由一系列等长/等行数的向量、因子、及/或矩阵组成，是一种特殊的list。数据框通常具有names和row.names两个附加属性，即为列名和行名。数据库是R里常用的数据结构类型。 可以用data.frame函数来创建，大的表格通常从硬盘中的表格文件读取。
(x <- data.frame(A=1:10, B=letters[1:10])) ## A B ## 1 1 a ## 2 2 b ## 3 3 c ## 4 4 d ## 5 5 e ## 6 6 f ## 7 7 g ## 8 8 h ## 9 9 i ## 10 10 j str(x) ## 'data.frame': 10 obs. of 2 variables: ## $ A: int 1 2 3 4 5 6 7 8 9 10 ## $ B: chr "a" "b" "c" "d" ... mode(x) ## [1] "list" class(x) ## [1] "data.frame" attributes(x) ## $names ## [1] "A" "B" ## ## $class ## [1] "data.frame" ## ## $row.names ## [1] 1 2 3 4 5 6 7 8 9 10 is.list(x) ## [1] TRUE is.data.frame(x) ## [1] TRUE (y <- unclass(x)) ## $A ## [1] 1 2 3 4 5 6 7 8 9 10 ## ## $B ## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" ## ## attr(,"row.names") ## [1] 1 2 3 4 5 6 7 8 9 10 str(y) ## List of 2 ## $ A: int [1:10] 1 2 3 4 5 6 7 8 9 10 ## $ B: chr [1:10] "a" "b" "c" "d" ... ## - attr(*, "row.names")= int [1:10] 1 2 3 4 5 6 7 8 9 10 mode(y) ## [1] "list" attributes(y) ## $names ## [1] "A" "B" ## ## $row.names ## [1] 1 2 3 4 5 6 7 8 9 10
函数
将一系列代码放在一起，就成为函数(function)，函数可以预设参数。将函数f的名称f作为命令得到的是f的内容，而要调用函数f必须要要用f()的形式，需要指定的参数放在括号里。
f <- function(){ print("A") print("Z") 1+2 } f ## function(){ ## print("A") ## print("Z") ## 1+2 ## } f() ## [1] "A" ## [1] "Z" ## [1] 3
函数里可以设置参数，参数可以在函数中的代码里引用：
kf <- function(X) sqrt(X) kf() ## Error in kf() : 缺少参数"X",也没有缺省值 kf(X=4) ## [1] 2
可以给定一个默认值，这样在调用函数的时候如果不给指定参数的值，函数就采用默认值。
kf <- function(X=4) sqrt(X) kf() ## [1] 2 kf(X=4) ## [1] 2 kf(X=9) ## [1] 3
可以设定很多的参数，在函数中采用：
f <- function(x="A", y="B", a=1, b=2){ print(x) print(y) a+b } f() ## [1] "A" ## [1] "B" ## [1] 3 f("a", "b", 3, 5) ## [1] "a" ## [1] "b" ## [1] 8 f(b=10, x="OK") ## [1] "OK" ## [1] "B" ## [1] 11
函数里的参数可以按顺序写，也可以不安顺序但是指明谁是谁。
R的函数分为向量化的和非向量化的，向量化的函数对参数中的向量逐个元素地进行运算，逐一返回运算结果组成新的向量。运算符也是一种函数，只是用法和一般的函数有别。在使用函数时需要区分向量化的和非向量化的区别。 一般来说，常用的一些函数用习惯了就知道怎么用了。
函数默认返回最后一个表达式的值，如果需要指定返回的值，使用return()函数。
f <- function(x=3, y=4){ X <- x^2 Y <- y^2 Y Z <- X + Y z <- Z^(1/2) } f()
当函数中没有return语句指定返回值的时候，如果最后一个命令也不是对象名时，那么函数将不会返回值。
f <- function(x=3, y=4){ X <- x^2 Y <- y^2 Z <- X + Y Z z <- Z^(1/2) z } f() ## [1] 5
当函数中没有return语句指定返回值的时候，如果最后一个命令是对象名时，那么函数将返此对象的值。
f <- function(x=3, y=4){ X <- x^2 Y <- y^2 Z <- X + Y return(Z) z <- Z^(1/2) z } f() ## [1] 25 f <- function(x=3, y=4){ X <- x^2 Y <- y^2 Z <- X + Y return(Z) z <- Z^(1/2) return(z) } f() ## [1] 25
当函数中有return语句指定返回值的时候，将返回第一个return返回的值。
对象类型 (II)
除了上边提到的比较基本的类型以外，还有一些在基本的类型上增加了不同的class和attributes的对象，这些对象有不同的用途。
因子
因子(factor) 是一个特殊的复合对象，可以看做特殊的向量，只能含有特定的值，用于存储分类数据。因子的模式是整数型的向量，具有有levels附加属性和factor class。 factor类型让因子与常规的整数型向量不同，而levels定义允许的值。用factor函数创建因子，将一个普通的向量转为因子,可以通过levels参数来改变levels及其顺序。
x <- c("家畜", "家畜", "家畜", "家禽", "家禽") (y <- factor(x)) ## [1] 家畜 家畜 家畜 家禽 家禽 ## Levels: 家畜 家禽 mode(y) ## [1] "numeric" class(y) ## [1] "factor" attributes(y) ## $levels ## [1] "家畜" "家禽" ## ## $class ## [1] "factor" ## str(y) ## Factor w/ 2 levels "家畜","家禽": 1 1 1 2 2 as.numeric(y) ## [1] 1 1 1 2 2 as.character(y) ## [1] "家畜" "家畜" "家畜" "家禽" "家禽"
y[6] <- "家禽" y[7] <- "野鸟" ## Warning message: ## In `[<-.factor`(`*tmp*`, 7, value = "野鸟") : 因子层次有错，产生了NA c(factor(c("a", "b")), factor(c("c", "d"))) ## [1] a b c d ## Levels: a b c d
不能增加levels里没有的元素，野鸟并不在y的levels之中，因此这里无法将y中的值换成野鸟，因而产生NA。可以用c将两个向量合并，levels也合并了
y <- factor(x, levels=c("家畜", "家禽")) order(y) ## [1] 1 2 3 4 5 sort(y) ## [1] 家畜 家畜 家畜 家禽 家禽 ## Levels: 家畜 家禽 z <- factor(x, levels=c("家禽", "家畜")) order(z) ## [1] 4 5 1 2 3 sort(z) ## [1] 家禽 家禽 家畜 家畜 家畜 ## Levels: 家禽 家畜 table(x) ## x ## 家畜 家禽 ## 3 2 table(factor(x, levels=c("家禽", "家畜"))) ## ## 家禽 家畜 ## 2 3 table(factor(x, levels=c("家禽", "家畜", "野鸟"))) ## ## 家禽 家畜 野鸟 ## 2 3 0
(z <- factor(c(1, 2, 4, 5, 7))) ## [1] 1 2 4 5 7 ## Levels: 1 2 4 5 7 as.numeric(z) ## [1] 1 2 3 4 5 as.character(z) ## [1] "1" "2" "4" "5" "7" as.numeric(as.character(z)) ## [1] 1 2 4 5 7 attr(z, "class") <- "numeric" as.character(z) ## [1] "1" "2" "3" "4" "5"
因子的本质是整数型的向量，只是因子的每个整数都对应着一个level，而level是字符型的。因此，如果把一个数值型向量转成因子，再由因子转回来的时候就要注意，可能得到的数字全变了。 如果强行将class改成numeric，这时候转换为character的时候就不按factor来了。
(X <- c(4, 6, 8)) ## [1] 4 6 8 str(X) ## num [1:3] 4 6 8 (Y <- as.factor(X)) ## [1] 4 6 8 ## Levels: 4 6 8 str(Y) ## Factor w/ 3 levels "4","6","8": 1 2 3 mode(Y) ## [1] "numeric" (Z <- as.numeric(Y)) ## [1] 1 2 3 str(Z) ## num [1:3] 1 2 3 (U <- as.character(Y)) ## [1] "4" "6" "8" str(U) ## chr [1:3] "4" "6" "8" (V <- as.numeric(U)) ## [1] 4 6 8 str(V) ## num [1:3] 4 6 8 (W <- as.numeric(as.character(Y))) ## [1] 4 6 8
在先前的一些版本的R中，默认情况下，在读入表格到data.frame或者创建data.frame的时候，会默认将字符型的向量转成因子，因此，当读数的数字一列中参杂字符的时候，原本应该当作数字读入的数据可能被转成了因子。 这时候如果要转成数值型的，就需要留意，以免出错。
日期与时间
Sys.time()得到当前的时间，得到的对象本质上还是数值型的，其class是“POSIXct” “POSIXt”。
(x <- Sys.time()) ## [1] "2022-12-05 00:54:00 CST" mode(x) ## [1] "numeric" class(x) ## [1] "POSIXct" "POSIXt" (x <- Sys.Date()) ## [1] "2022-12-05" mode(x) ## [1] "numeric" class(x) ## [1] "Date" as.Date("1949-10-01") ## [1] "1949-10-01" Sys.Date() - as.Date("1949-10-01") ## Time difference of 26728 days as.Date("2049-10-01") - Sys.Date() ## [1] "2022-12-10" as.Date("25/10/2022", format="%d/%d/%Y") (X <- Sys.Date()) ## [1] "2022-12-05" (Y <- as.numeric(x)) ## [1] 19331 X - Y ## [1] "1970-01-01" (Z <- as.Date(Y, origin="1970-01-01")) ## [1] "2022-12-05"
默认情况下转化为数字是从“1970-01-01”算第一天的。从数字转为日期的时候，一定要正确设置origin。从EXCEL读入日期数据的时候，得到的是数值，这时候就需要转换成日期，需要注意的是Windows 下EXCEL日期origin=“1899-12-30”， Mac的EXCEL的为“1904-01-01”，如果不明确文件来源的时候，转换之后需要核对一下是否正确。
table
table将统计因子中levels各个分类的数目。如果给定的是可以转化为因子的向量，那么table将之转化为因子并作出统计。
set.seed(2022) (x <- sample(LETTERS[1:5], 10, replace=TRUE)) ## [1] "D" "C" "C" "D" "D" "C" "A" "E" "B" "E" (y <- table(x)) ## x ## A B C D E ## 1 1 3 3 2 mode(y) ## [1] "numeric" class(y) ## [1] "table" str(y) ## 'table' int [1:5(1d)] 1 1 3 3 2 ## - attr(*, "dimnames")=List of 1 ## ..$ x: chr [1:5] "A" "B" "C" "D" ... y[x] ## x ## D C C D D C A E B E ## 3 3 3 3 3 3 1 2 1 2
对象部分提取或替换
前边已经说到了对象取子集，取出对象的部分元素可以通过[、[[、$来，同时也可以对这部分进行赋值从而将其替换。基本用法如下：
x[i] x[i, j, ... , drop = TRUE] x[[i, exact = TRUE]] x[[i, j, ..., exact = TRUE]] x$name getElement(object, name)
取出对象x的部分，取子集不改变x，其中，向量i、j、……代表行、列、……（第三维、第四维、……），其取值可以是如下三种值：
1. 使用位置向量，向量前边加“-“表示不包括
2. 使用等长的逻辑向量，TRUE对应需要取出的
3. 使用元素名称向量，当元素无名称时就不好用了
x[i] <- value x[i, j, ...] <- value x[[i]] <- value x$name <- value
可以通过赋值改变相应的部分，这时候会改变x。
用如下命令了解更多关于对象部分提取或替换 (Extract or Replace Parts of an Object)的用法：
? "["
向量取子集
使用索引序号(index)来取出向量的部分，依次按照向量中的序号取出：
X <- c(1, 3, 5, 7, 9, 2, 4, 6,8, 10) names(X) <- LETTERS[1:10] X ## A B C D E F G H I J ## 1 3 5 7 9 2 4 6 8 10 X[] ## A B C D E F G H I J ## 1 3 5 7 9 2 4 6 8 10 X[0] ## named numeric(0) X[1:4] ## A B C D ## 1 3 5 7 X[c(1, 2, 4, 3, 3)] ## A B D C C ## 1 3 7 5 5 X[-c(1:4)] ## E F G H I J ## 9 2 4 6 8 10 X[-1:4] ## Error in X[-1:4] : 只有负下标里才能有零 X[c(1,3,5,8,6,5,4,2,50)] ## A C E H F E D B <NA> ## 1 5 9 6 2 9 7 3 NA X[-c(1,3,5,8,6,5,4,2,50)] ## G I J ## 4 8 10
保持空的时候得到的就是对象本身。index是从1开始的，给0就得到长度为0向量。
X[1, 2, 4] ## Error in X[1, 2, 4] : 量度数目不对
这里的索引需要给一个向量，如上这样就是给三个向量了，对于向量取子集来说就会报错了。这是比较常见的一种错误，需要注意。
根据索引序号来取子集可以用于对向量进行排序：
X[order(X)] ## A F B G C H D I E J ## 1 2 3 4 5 6 7 8 9 10 X[order(X, decreasing=TRUE)] ## J E I D H C G B F A ## 10 9 8 7 6 5 4 3 2 1
order这个函数得到的是对给定向量进行排序时取子集的序号的向量。
Y <- c(3, 5, 1, 7) order(Y) ## [1] 3 1 2 4 order(Y, decreasing = TRUE) ## [1] 4 2 1 3
对Y进行升序排序时，放在第一位的是原来的第3为，放在第二位的是原来的第2位，依次类推。默认是升序排列，降序排序需要加上参数 decreasing = TRUE
使用逻辑向量来取，对应TRUE的位置的元素被取出：
X > 5 ## A B C D E F G H I J ## FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE X[X > 5] ## D E H I J ## 7 9 6 8 10
使用元素名称向量取子集，对应名称被取出：
X[c("B", "F", "J")] ## B F J ## 3 2 10
这样使得向量可以被当做字典来用。
取子集的时候赋值，可以改变向量中的部分元素：
X <- 1:10 names(X) <- LETTERS[1:10] X ## A B C D E F G H I J ## 1 2 3 4 5 6 7 8 9 10 X[5] <- 1000 X ## A B C D E F G H I J ## 1 2 3 4 1000 6 7 8 9 10 X[1:5] <- 11:15 X ## A B C D E F G H I J ## 11 12 13 14 15 6 7 8 9 10 X[c("B", "C")] <- c("b", "c", "d") ## Warning message: ## In X[c("B", "C")] <- c("b", "c", "d") : 被替换的项目不是替换值长度的倍数 X ## A B C D E F G H I J ## "11" "b" "c" "14" "15" "6" "7" "8" "9" "10"
矩阵取子集
矩阵是二维的，可以通过行列的序号来取：
M <- matrix(1:20, ncol=5, byrow=TRUE, dimnames=list(letters[1:4], LETTERS[1:5])) M ## A B C D E ## a 1 2 3 4 5 ## b 6 7 8 9 10 ## c 11 12 13 14 15 ## d 16 17 18 19 20 M[1:2, c(1, 3, 5)] ## A C E ## a 1 3 5 ## b 6 8 10 M[4, 4] ## [1] 19
可以通过行列的名称来取:
M[c("a", "c"), c("D", "E")] ## D E ## a 4 5 ## c 14 15
可以通过行列等长的逻辑向量来取：
M[c(TRUE, FALSE, TRUE, FALSE), c(FALSE, TRUE, TRUE, FALSE, FALSE)] ## B C ## a 2 3 ## c 12 13 M[rownames(M) %in% c("a", "B", "c"), colnames(M) %in% c("a", "B", "c", "D")] ## B D ## a 2 4 ## c 12 14
同样可以通过赋值来对其部分进行替换：
M[1:2, c(1, 3, 5)] <- 0 M ## A B C D E ## a 0 2 0 4 0 ## b 0 7 0 9 0 ## c 11 12 13 14 15 ## d 16 17 18 19 20 M[1:2, 1:2] <- matrix(101:104, ncol=2, byrow=TRUE) M ## A B C D E ## a 101 102 0 4 0 ## b 103 104 0 9 0 ## c 11 12 13 14 15 ## d 16 17 18 19 20 M[3:4, 3:4] <- NA M ## A B C D E ## a 101 102 0 4 0 ## b 103 104 0 9 0 ## c 11 12 NA NA 15 ## d 16 17 NA NA 20
矩阵也可以本质上是向量，可以用一维的方式来进行取值和替换，这时候序号是像第一行然后第二行这样依次排序：
M[is.na(M)] <- 50 M ## A B C D E ## a 101 102 0 4 0 ## b 103 104 0 9 0 ## c 11 12 50 50 15 ## d 16 17 50 50 20 which(M > 50) [1] 1 2 5 6 11 12 15 16
比较实用的情况就是对NA值进行替换。
列表取子集
列表当做一维的来看，可以用方括号[]取子集，也可以用$或者[[]]取出其中的元素。
Y <- list(A=1:3, B=letters[1:5], C=LETTERS[6:10]) Y ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e" ## ## $C ## [1] "F" "G" "H" "I" "J" Y[1] ## $A ## [1] 1 2 3 Y[c("A","B")] ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e" Y[1:2] ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e"
可以使用$加元素的名称来取出列表的元素，$后边时元素的名称:
Y$A ## [1] 1 2 3 Y$B ## [1] "a" "b" "c" "d" "e" Y$A[2] ## [1] 2 Y$B[2:3] ## [1] "b" "c"
双重方括号[[]]可用于取出单一的元素:
Y[[1]] ## [1] 1 2 3 Y[["B"]] ## [1] "a" "b" "c" "d" "e"
请注意Y[1]和Y[[1]]的区别，Y["B"]和Y[["B"]]。
[[]]取元素可递归：
Z <- list( a = list( b = 9, c = "hello", d = list( e = "3rd", f = list( g="4th"))), h = 1:5) Z ## $a ## $a$b ## [1] 9 ## ## $a$c ## [1] "hello" ## ## $a$d ## $a$d$e ## [1] "3rd" ## ## $a$d$f ## $a$d$f$g ## [1] "4th" Z[[c("a", "b")]] ## [1] 9 Z[[c("a", "b", "c")]] ## Error in Z[[c("a", "b", "c")]] : 下标出界 Z[[c("a", "d", "f")]] ## $g ## [1] "4th" Z[[c("a", "d", "f", "g")]] ## [1] "4th" Z$a$d$f$g ## [1] "4th" Y[[c("A","B")]] ## Error in Y[[c("A", "B")]] : 下标出界 Y[c(1, 2)] ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e" Y[[c(1, 2)]] ## [1] 2
Z[[c("a", "b")]]递归就是先取Z的元素a，再取a的元素b，依次往后推。Z[[c("a", "d", "f")]]中选取出a，然后取出a中的d，然后d中没有f，所以报错说下标出界。其实这种递归的方式并没有连用$的方式直观和清晰。要留意Y[c(1, 2)]和Y[[c(1, 2)]]的区别。
Y <- list(A=1:3, B=letters[1:5], C=LETTERS[6:10]) Y ## $A ## [1] 1 2 3 ## ## $B ## [1] "a" "b" "c" "d" "e" ## ## $C ## [1] "F" "G" "H" "I" "J" Y[1:2] <- list(D=1:5, E=10:6) Y
数据框取子集
数据框本质上是列表，适合于列表的取子集的方式也适合于数据框，同时，数据框又像矩阵一样是二维的，可以通过行和列来取：
X <- data.frame( N = 11:16, S = c("a", "b", "a", "c", "a", "b")) row.names(X) <- LETTERS[1:6] X ## N S ## A 11 a ## B 12 b ## C 13 a ## D 14 c ## E 15 a ## F 16 b
可以通过$或者[[]]来取出某列：
X$N ## [1] 11 12 13 14 15 16 X[["S"]] ## [1] "a" "b" "a" "c" "a" "b"
可以通过行列的序号来取出部分：
X[c(1, 3, 6), ] ## N S ## A 11 a ## C 13 a ## F 16 b
利用序号取子集，可用于对数据框进行排序：
X[order(X$N, decreasing=TRUE), ] ## N S ## F 16 b ## E 15 a ## D 14 c ## C 13 a ## B 12 b ## A 11 a X[order(factor(X$S, levels=c("b", "a", "c"))), ] ## N S ## B 12 b ## F 16 b ## A 11 a ## C 13 a ## E 15 a ## D 14 c X <- X[order(factor(X$S, levels=c("b", "a", "c"))), ] X ## N S ## B 12 b ## F 16 b ## A 11 a ## C 13 a ## E 15 a ## D 14 c
通过逻辑向量来取：
X$N < 15 & X$N > 11 ## [1] FALSE TRUE TRUE TRUE FALSE FALSE X[c(FALSE, TRUE, TRUE, TRUE, FALSE, FALSE), ] ## N S ## B 12 b ## C 13 a ## D 14 c X[X$N < 15 & X$N > 11, ] ## N S ## B 12 b ## C 13 a ## D 14 c X[X$S %in% c("b", "c"), ] ## N S ## B 12 b ## D 14 c ## F 16 b X[X$S == "a", ] ## N S ## A 11 a ## C 13 a ## E 15 a X[ , names(X) %in% c("S")] ## [1] "b" "b" "a" "a" "a" "c" X[X$N < 5, names(X) %in% c("S")] ## character(0)
通过行列的名称来取子集：
X[c("D", "E", "A"), ] ## N S ## D 14 c ## E 15 a ## A 11 a X[ , c("S", "N")] ## S N ## B b 12 ## F b 16 ## A a 11 ## C a 13 ## E a 15 ## D c 14 X[ , c(2, 1)] ## S N ## B b 12 ## F b 16 ## A a 11 ## C a 13 ## E a 15 ## D c 14
当取单行或者单列的时候，需要留意取出来的时候是不是对象类型改变了。用[]取单列的时候，默认得到的是本列的向量，而不是一个单列的数据框。
X[ , "N"] ## [1] 11 12 13 14 15 16 X[ , "N", drop=TRUE] ## [1] 11 12 13 14 15 16 X[ , "N", drop=FALSE] ## N ## A 11 ## B 12 ## C 13 ## D 14 ## E 15 ## F 16 X["A", ] ## N S ## A 11 a X["A", , drop=FALSE] ## N S ## A 11 a X["A", , drop=TRUE] ## $N ## [1] 11 ## ## $S ## [1] "a"
取子集与赋值合用，实现对象中子集值的改变：
X <- data.frame(N=1:6, S=c("a", "b", "a", "c", "a", "b")) X ## N S ## 1 1 a ## 2 2 b ## 3 3 a ## 4 4 c ## 5 5 a ## 6 6 b X[X$S == "a", ]$N <- 10 X ## N S ## 1 10 a ## 2 2 b ## 3 10 a ## 4 4 c ## 5 10 a ## 6 6 b X$X <- "Y" X ## N S X ## 1 10 a Y ## 2 2 b Y ## 3 10 a Y ## 4 4 c Y ## 5 10 a Y ## 6 6 b Y X$Y <- 1:10 ## Error in `$<-.data.frame`(`*tmp*`, Y, value = 1:10) : ## 替换数据里有10行，但数据有6 X$Y <- 11:16 X ## N S X Y ## 1 10 a Y 11 ## 2 2 b Y 12 ## 3 10 a Y 13 ## 4 4 c Y 14 ## 5 10 a Y 15 ## 6 6 b Y 16
同样，有时候要对整个数据框里的NA值进行替换，也可以像矩阵里替换一样：
X <- data.frame(A=1:6, B=letters[1:6], C=LETTERS[1:6], D=6:1) X ## A B C D ## 1 1 a A 6 ## 2 2 b B 5 ## 3 3 c C 4 ## 4 4 d D 3 ## 5 5 e E 2 ## 6 6 f F 1 X$A[c(1, 3, 5)] <- NA X$C[c(2, 4, 6)] <- NA X ## A B C D ## 1 NA a A 6 ## 2 2 b <NA> 5 ## 3 NA c C 4 ## 4 4 d <NA> 3 ## 5 NA e E 2 ## 6 6 f <NA> 1 which(is.na(X)) ## [1] 1 3 5 14 16 18 X[is.na(X)] <- 0 X ## A B C D ## 1 0 a A 6 ## 2 2 b 0 5 ## 3 0 c C 4 ## 4 4 d 0 3 ## 5 0 e E 2 ## 6 6 f 0 1
需要注意，赋值时，不能把sample函数直接放到Index的位置，这样得不到预期：
X <- data.frame(N=1:6, S=c("a", "b", "a", "c", "a", "b"), U=LETTERS[1:6]) X ## N S U ## 1 1 a A ## 2 2 b B ## 3 3 a C ## 4 4 c D ## 5 5 a E ## 6 6 b F set.seed(1); sample(1:6, 3) ## [1] 1 4 3 set.seed(1); X[sample(1:6, 3), ]$N <- 8:10 X ## N S U ## 1 8 a A ## 2 9 c D ## 3 3 a C ## 4 4 c D ## 5 5 a E ## 6 10 a C set.seed(1); I <- sample(1:6, 3) X[I, ]$N <- 8:10 X ## N S U ## 1 8 a A ## 2 9 c D ## 3 10 a C ## 4 9 c D ## 5 5 a E ## 6 10 a C
如果是要随机改变某列的一些行，那么需要先把随机数产生出来，再放到Index的位置。
data.frame是一个特殊的列表，通过如下例子，也能够更好地理解[]和[[]]的区别：
D <- data.frame(X=1:2, Y=3:4) D ## X Y ## 1 1 3 ## 2 2 4 D["X"] ## X ## 1 1 ## 2 2 D[c("X", "Y")] ## X Y ## 1 1 3 ## 2 2 4 D$X ## [1] 1 2 D[["X"]] ## [1] 1 2 class(D["X"]) ## [1] "data.frame" class(D[["X"]]) ## [1] "integer"
可以看到，D[“X”]得到的依然还是一个数据框，含有X一列，而D[["X"]]得到的是X这一列，是D的元素，是个向量。
数据库取子集(应用示例)
CJ <- read.delim("ftp://ftp.biolab.wang/data/各科成绩.tsv", fileEncoding="GBK") head(CJ) ## 学号 姓名 专业 班级 性别 科目 成绩 ## 1 101200101 罗诗超 专业01 2001班 男 细胞生物学 57 ## 2 101200102 韦海凤 专业01 2001班 女 细胞生物学 77 ## 3 101200103 徐钧月 专业01 2001班 男 细胞生物学 71 ## 4 101200104 黄妮亨 专业01 2001班 女 细胞生物学 85 ## 5 101200105 欧炳彤 专业01 2001班 男 细胞生物学 64 ## 6 101200106 潘叶 专业01 2001班 女 细胞生物学 73 str(CJ) ## 'data.frame': 4848 obs. of 7 variables: ## $ 学号: int 101200101 101200102 101200103 101200104 101200105 101200106 101200107 101200108 101200109 101200110 ... ## $ 姓名: chr "罗诗超" "韦海凤" "徐钧月" "黄妮亨" ... ## $ 专业: chr "专业01" "专业01" "专业01" "专业01" ... ## $ 班级: chr "2001班" "2001班" "2001班" "2001班" ... ## $ 性别: chr "男" "女" "男" "女" ... ## $ 科目: chr "细胞生物学" "细胞生物学" "细胞生物学" "细胞生物学" ... ## $ 成绩: int 57 77 71 85 64 73 77 78 81 56 ...
拿到一个数据框，首先通过head和str看看基本的信息。
CJ2001 <- CJ[CJ$班级 == "2001班", ] head(CJ2001) ## 学号 姓名 专业 班级 性别 科目 成绩 ## 1 101200101 罗诗超 专业01 2001班 男 细胞生物学 57 ## 2 101200102 韦海凤 专业01 2001班 女 细胞生物学 77 ## 3 101200103 徐钧月 专业01 2001班 男 细胞生物学 71 ## 4 101200104 黄妮亨 专业01 2001班 女 细胞生物学 85 ## 5 101200105 欧炳彤 专业01 2001班 男 细胞生物学 64 ## 6 101200106 潘叶 专业01 2001班 女 细胞生物学 73 CJ2001 <- CJ2001[order(CJ2001$成绩, decreasing = TRUE), ] head(CJ2001) ## 学号 姓名 专业 班级 性别 科目 成绩 ## 1238 101200126 刘丽 专业01 2001班 女 体育 94 ## 1844 101200126 刘丽 专业01 2001班 女 微生物学 94 ## 4268 101200126 刘丽 专业01 2001班 女 分子生物学 94 ## 1822 101200104 黄妮亨 专业01 2001班 女 微生物学 93 ## 1225 101200113 黄建华 专业01 2001班 女 体育 91 ## 1216 101200104 黄妮亨 专业01 2001班 女 体育 90 tail(CJ2001) ## 学号 姓名 专业 班级 性别 科目 成绩 ## 2480 101200156 江雨宇 专业01 2001班 男 生物化学 46 ## 661 101200155 黄佳涛 专业01 2001班 女 解剖学 45 ## 649 101200143 桂泳 专业01 2001班 女 解剖学 43 ## 1874 101200156 江雨宇 专业01 2001班 男 微生物学 43 ## 4298 101200156 江雨宇 专业01 2001班 男 分子生物学 43 ## 3691 101200155 黄佳涛 专业01 2001班 女 生理学 42 CJ[CJ$姓名 == "陆双全", ] ## 学号 姓名 专业 班级 性别 科目 成绩 ## 47 101200147 陆双全 专业01 2001班 女 细胞生物学 61 ## 653 101200147 陆双全 专业01 2001班 女 解剖学 60 ## 1259 101200147 陆双全 专业01 2001班 女 体育 68 ## 1865 101200147 陆双全 专业01 2001班 女 微生物学 64 ## 2471 101200147 陆双全 专业01 2001班 女 生物化学 64 ## 3077 101200147 陆双全 专业01 2001班 女 遗传学 63 ## 3683 101200147 陆双全 专业01 2001班 女 生理学 66 ## 4289 101200147 陆双全 专业01 2001班 女 分子生物学 57 CJ[CJ$姓名 == "陆双全" & CJ$成绩 < 60, ]$成绩 <- 60 CJ[CJ$姓名 == "陆双全", ] ## 学号 姓名 专业 班级 性别 科目 成绩 ## 47 101200147 陆双全 专业01 2001班 女 细胞生物学 61 ## 653 101200147 陆双全 专业01 2001班 女 解剖学 60 ## 1259 101200147 陆双全 专业01 2001班 女 体育 68 ## 1865 101200147 陆双全 专业01 2001班 女 微生物学 64 ## 2471 101200147 陆双全 专业01 2001班 女 生物化学 64 ## 3077 101200147 陆双全 专业01 2001班 女 遗传学 63 ## 3683 101200147 陆双全 专业01 2001班 女 生理学 66 ## 4289 101200147 陆双全 专业01 2001班 女 分子生物学 60
R的运算
向量运算符
算术运算：加(+)、减(-)、乘(*)、除以(/)、指数(^)、整除(%/%)、取模(%%)
比较运算：相等(==)、不等(!=)、大于(>)、小于(<)、大于等于(>=)、小于等于(<=)
逻辑运算：或(|)、且(&)、非(!)
算术运算针对数值型或者复数型的对象进行加减乘数的基本算术，比较运算、逻辑运算结果为逻辑型。如：
10 + 2*3 + 4/5 -2*6 + 2^3 ## [1] 12.8 11.3/2 ## [1] 5.65 11.3 %/% 2 ## [1] 5 11.3 %% 2 ## [1] 1.3 1 == 2 ## [1] FALSE 1 <= 2 ## [1] TRUE 1 <= 2 & 1 == 2 ## [1] FALSE 1 <= 2 | 1 == 2 ## [1] TRUE !(1 <= 2 & 1 == 2) ## [1] TRUE
向量运算
如上提到的运算符都可以用于向量和向量之间，两个向量之前相同位置上的数值之间进行运算，运算结果按顺序返回生成新的向量。当两个向量长度不等时，较短的向量将重新从第一个元素开始重复使用，如此直到较长的向量所有元素都完成。这样当一个向量与长度为1的向量进行运算时，相当于每个元素都与这个值做同样的运算。
1:5 + 10:14 ## [1] 11 13 15 17 19 1:5 + 10:19 ## [1] 11 13 15 17 19 16 18 20 22 24 1:5 + 10:23 ## [1] 11 13 15 17 19 16 18 20 22 24 21 23 25 27 ## Warning message: ## In 1:5 + 10:23 : 长的对象长度不是短的对象长度的整倍数
更多运算符
&& (且，非向量化的，用于长度为1的向量)， || (或，非向量化的，用于长度为1的向量) ，一般用于 if() 中。
%in% (前一个向量的元素是否出现在后一个向量里,对前一个向量中每个元素依次判断返回向量结果)
<- (向左赋值)
-> (向右赋值)
TRUE && c(TRUE, FALSE) ## [1] TRUE ## Warning message: ## In TRUE && c(TRUE, FALSE) : 'length(x) = 2 > 1' in coercion to 'logical(1)' TRUE && c(FALSE, TRUE) ## [1] FALSE ## Warning message: ## In TRUE && c(FALSE, TRUE) : 'length(x) = 2 > 1' in coercion to 'logical(1)' FALSE || c(TRUE, FALSE) ## [1] TRUE ## Warning message: ## In FALSE || c(TRUE, FALSE) : ## 'length(x) = 2 > 1' in coercion to 'logical(1)' FALSE || c(FALSE, TRUE) ## [1] FALSE ## Warning message: ## In FALSE || c(FALSE, TRUE) : ## 'length(x) = 2 > 1' in coercion to 'logical(1)' c(1,3,5) %in% c(2,3,5,6) ## [1] FALSE TRUE TRUE c(2,3,5,6) %in% c(1,3,5) ## [1] FALSE TRUE TRUE FALSE y <- 1 -> x y ## [1] 1 x ## [1] 1
管道操作符
我个人是不太建议使用%>%或|>这样的管道的，我认为还是经典的R的向左赋值更清晰明了。然而，管道却被不少人吹捧，因此，这里还是提一下，以便大家见到了能看懂。
管道操作符将其前边的表达式得到的值传递给它后边的表达式或者函数，作为后边函数的第一个参数，或者将值传到占位符（%>%对应.而|>对应_）所在的位置。
%>%
magrittr包使用管道操作符%>%等。使用install.packages('magrittr')命令安装，library(magrittr)加载后使用：
x %>% f(y) 等同于 f(x, y)，而x %>% f(y, arg=.) 等同于 f(y, arg=x)
1:10 %>% mean() ## [1] 5.5 1:5 %>% round(3.1415926, digits=.) ## [1] 3.10000 3.14000 3.14200 3.14160 3.14159 1:10 %>% sd ## [1] 3.02765 CJ <- read.delim("ftp://ftp.biolab.wang/data/各科成绩.tsv", fileEncoding="GBK") CJ %>% subset(姓名=="罗诗超" & 科目=="细胞生物学") ## 学号 姓名 专业 班级 性别 科目 成绩 ## 1 101200101 罗诗超 专业01 2001班 男 细胞生物学 57
|>
现在R的base包已有其管道操作符|>, 用法与%>%一样，不一样的地方是|>后边不能只给函数名，而是需要是函数调用，也就是说括号不能省略。另外，占位符用的是下划线_。
1:10 |> sqrt() ## [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427 3.000000 3.162278 1:5 |> round(3.01234567890, digits=_) ## [1] 3.00000 3.01000 3.01200 3.01230 3.01235
流程控制
流程控制结构控制代码的执行流程，主要包括循环和条件判断。循环用于重复执行一段代码，直到满足某个条件才终止。在一个循环中，对对重复执行的这一段代码的一次执行叫做一次迭代。通过多次迭代重复执行某个过程或者步骤，逐步改进或者逼近目标。常用的就是流程控制为for, repeat, while等，条件判断有if, if … else …等，终止一个循环用break，终止一次迭代用next。
for 循环
对向量中的每个元素依次执行操作，执行的操作可以不用到这些元素。比如：
for(i in c(0, 11, 235, 812, 2134)){ print(i) } ## [1] 0 ## [1] 11 ## [1] 235 ## [1] 812 ## [1] 2134 for(i in c(0, 11, 235, 812, 2134)){ I <- sprintf("%04d", i) print(I) } ## [1] "0000" ## [1] "0011" ## [1] "0235" ## [1] "0812" ## [1] "2134"
如上依次进行格式化之后打印，当然还能有别的很多你自己希望的步骤，这里的i就是依次取用向量中的值。
for(i in c(0, 11, 235, 812, 2134)){ I <- sample(base::LETTERS, 1) print(I) } ## [1] "B" ## [1] "A" ## [1] "G" ## [1] "J" ## [1] "B"
如上，操作可以不用到向量中的元素，这时候其实就是个简单的计数，运行向量长度那么多次。LETTERS是base包中内置的向量，元素为A-Z这26个大写的字母。sample用于从中随机抽出一个，大家练习的时候可能得到的不是这些字母。
在R中的for循环速度比较慢，有时候需要避免，而是采用向量化的运算比较合适。这个用system.time比较一下可以看出来:
# 如下是向量化的运算，可以看到，速度是比较快的 system.time({I <- sum((1:100000000)^2); print(I)}) # 如下用for来算，很显然使用了几十上百倍的时间 system.time({I <- 0; for(i in c(1:100000000)) I <- I + i^2; print(I)})
同样是对1~100000000这一亿个数进行求平方，用for来一个个算就比向量化的计算要慢得多。
while 循环
while(当)满足条件时执行，比如
X <- c(0, 11, 235, 812, 2134) i <- 1 while(i <= 4){ print(X[i]) i <- i+1 } ## [1] 0 ## [1] 11 ## [1] 235 ## [1] 812
如上依次打印到屏幕，直到i大于5就停止了。 第一次判断，i是1，比4小，于是运行第一次迭代，运行之后，i变成了2，依然满足条件，于是第二次迭代，当进行了4次迭代之后，i变成了5，不再满足条件，就退出循环不再迭代了，所以X的第5个元素并没有被print出来。
X <- c(0, 11, 235, 812, 2134) i <- 1 while(i <= 4){ I <- sprintf("%04d", X[i]) print(I) i <- i+1 } ## [1] "0000" ## [1] "0011" ## [1] "0235" ## [1] "0812"
如上依次进行格式化之后打印，当然还能有别的很多你自己希望的步骤，这里的i就是依次取用向量中的值。
很显然，这里while也能实现想要用for来做的事情，但不一定比for做得好。比如，我们如果要把1~10000000000这一百亿个数的平方和，那么试试如下三种做法，体会一下:
system.time({S <- sum((1:10000000000)^2); print(S)}) # 如上方法错误，比如我的电脑上显示错误为：“错误: 无法分配大小为74.5 Gb的矢量” system.time({S <- 0; for(i in 1:10000000000){S <- S + i^2}; print(S)}) # 依然能够运行 system.time({S <-0; i <- 1; while(i <= 10000000000){S <- S + i^2; i <- i + 1}; print(S)}) # 依然能够运行，但是时间比如for要长一些，大概是因为这里每次迭代都多了一步i <- i+1的缘故
repeat 循环
进行无条件的循环，因此需要在流程中制定退出机制。使用if和break决定什么时候退出：
X <- c(0, 11, 235, 812, 2134) i <- 1 repeat{ print(X[i]) if(i >= 5) break i <- i+1 } ## [1] 0 ## [1] 11 ## [1] 235 ## [1] 812 ## [1] 2134
每次repeat之后i的值都是变化的。注意break在不同的位置结果可能是不同的，需要留意程序中各个值的变化：
X <- c(0, 11, 235, 812, 2134) i <- 1 repeat{ print(X[i]) i <- i+1 if(i >= 5) break } ## [1] 0 ## [1] 11 ## [1] 235 ## [1] 812
if 与 if...else...
当满足条件时执行一次, 后边可以用(也可以不用)else来指明不满足条件时执行的操作。 如下：
if(1 < 2 ) print("TRUE") ## [1] "TRUE" if(1 > 2 ) print("TRUE") if(1 < 2 ) print("TRUE") else print("Not TRUE") ## [1] "TRUE" if(1 > 2 ) print("TRUE") else print("Not TRUE") ## [1] "Not TRUE"
else 后边可以紧接着用if，这样就可以有很多的判断，很多的分支。如下这样，把TRUE换成适当的判断语句，把…换成需要执行的命名：
if(TRUE){ ... }else if(TRUE){ ... }else if(TRUE){ ... }else if(TRUE){ ... }esle{ ... }
比如：
i=1 while(i <= 20){ if(i < 5){ print(sprintf("i=%s, i < 5", i)) i = i + 2 } else if(i < 8) { print(sprintf("i=%s, i >= 5 & i < 8", i)) i = i + 1 } else { print(sprintf("i=%s, i >= 8", i)) i = i + 3 } } ## [1] "i=1, i < 5" ## [1] "i=3, i < 5" ## [1] "i=5, i >= 5 & i < 8" ## [1] "i=6, i >= 5 & i < 8" ## [1] "i=7, i >= 5 & i < 8" ## [1] "i=8, i >= 8" ## [1] "i=11, i >= 8" ## [1] "i=14, i >= 8" ## [1] "i=17, i >= 8" ## [1] "i=20, i >= 8"
满足条件的时候执行，需预先定义i，且循环中i应变化，否则即是死循环。
for(i in 1:10){ if(i < 5){ print(sprintf("i=%s, i < 5", i)) } else if(i < 8) { print(sprintf("i=%s, i >= 5 & i < 8", i)) } else { print(sprintf("i=%s, i >= 8", i)) } } ## [1] "i=1, i < 5" ## [1] "i=2, i < 5" ## [1] "i=3, i < 5" ## [1] "i=4, i < 5" ## [1] "i=5, i >= 5 & i < 8" ## [1] "i=6, i >= 5 & i < 8" ## [1] "i=7, i >= 5 & i < 8" ## [1] "i=8, i >= 8" ## [1] "i=9, i >= 8" ## [1] "i=10, i >= 8"
for对所有的元素依次执行，if做出判断，else if 部分可以没有，也可以有多个。
for(i in 1:10){ if(i < 5){ next print(sprintf("%s < 5", i)) } else if(i < 8) { print(sprintf("%s >= 5 & %s < 8", i, i)) } else { print(sprintf("%s >= 8", i)) break } } ## [1] "5 >= 5 & 5 < 8" ## [1] "6 >= 5 & 6 < 8" ## [1] "7 >= 5 & 7 < 8" ## [1] "8 >= 8"
next是跳出循环的当前迭代，进入下一次迭代，而break是终止整个当前的循环，这里是for循环。
ifelse
根据一个逻辑向量中的TRUE和FALSE对两个向量中的值进行选择，同样返回一个向量。比如：
ifelse(1:10 <= 5, "小", "大") ## [1] "小" "小" "小" "小" "小" "大" "大" "大" "大" "大" ifelse(c(1, 3, 5, 7, 9, 2, 4, 6, 8, 10) <= 5, "A", "B") ## [1] "A" "A" "A" "B" "B" "A" "A" "B" "B" "B" ifelse(1:10 <= 5, base::LETTERS, base::letters) ## [1] "A" "B" "C" "D" "E" "f" "g" "h" "i" "j" ifelse(1:10 <= 5, base::LETTERS[1:2], base::letters[1:3]) ## [1] "A" "B" "A" "B" "A" "c" "a" "b" "c" "a"
给ifelse的三个参数都是向量，如果后边的两个向量长度小于第一个向量(逻辑向量)的长度，则重复使用，直到与第一个向量等长。
ifelse(c(1, 3, 5, 7, 9, 2, 4, 6, 8, 10) <= 5, base::LETTERS[1:2], base::letters[1:3]) ## [1] "A" "B" "A" "a" "b" "B" "A" "b" "c" "a"
比如上边这个例子，将后边的两个向量重复之后与第一个向量的值对应如下,这样就很好理解结果了：
TRUE A a TRUE B b TRUE A c FALSE B a FALSE A b TRUE B c TRUE A a FALSE B b FALSE A c FALSE B a
switch
根据表达式(表达式的运算结果是一个数字或者字符串)来选择返回的值，使用? switch查看更多的。比如下边这些例子：
for(X in c("A", "CC", "DD")){ Y <- switch(EXPR=substr(X, 1, 1), A="a", B="bb", C="ccc") print(Y) } ## [1] "a" ## [1] "ccc" ## NULL
switch用于根据EXPR的值来选择，当这个值并没有对应的值的时候，得到NULL。
for(X in c("A", "C", "D", "E")){ Y <- switch(EXPR=X, A="a", B="bb", C="ccc", "dddd") print(Y) } ## [1] "a" ## [1] "ccc" ## [1] "dddd" ## [1] "dddd"
而如果存在并没有名称的选项的时候，所有的没有对应项的值都会得到这个没有名称的值。
for(X in 1:5){ Y <- switch(EXPR=X, A="a", B="bb", C="ccc", "dddd") print(Y) } ## [1] "a" ## [1] "bb" ## [1] "ccc" ## [1] "dddd" ## NULL
可以使用数字，当使用数字的时候，只是和顺序有关，而和选项是否有名称并没有多大的关系。
for(X in 0:10){ Y <- switch(EXPR=sqrt(X), A="a", "dddd") print(Y) } ## NULL ## [1] "a" ## [1] "a" ## [1] "a" ## [1] "dddd" ## [1] "dddd" ## [1] "dddd" ## [1] "dddd" ## [1] "dddd" ## NULL ## NULL for(X in 0:10){ Y <- switch(EXPR=sqrt(X), A="a", "dddd") print(sprintf("EXPR=%f, Y gets %s", sqrt(X), Y)) } ## character(0) ## [1] "EXPR=1.000000, Y gets a" ## [1] "EXPR=1.414214, Y gets a" ## [1] "EXPR=1.732051, Y gets a" ## [1] "EXPR=2.000000, Y gets dddd" ## [1] "EXPR=2.236068, Y gets dddd" ## [1] "EXPR=2.449490, Y gets dddd" ## [1] "EXPR=2.645751, Y gets dddd" ## [1] "EXPR=2.828427, Y gets dddd" ## character(0) ## character(0)
可以看到，当EXPR的值不是整数的时候，按照整数部分来取值。这些也能够用if…else…来实现。
break 与 next
break用于退出for、while或repeat循环，不再进行后续的迭代，next用于结束当前的迭代，进入下一次迭代。当多个循环套用时，break和next都只作用于所在的循环。使用时要注意终止的位置。
for(i in 1:10){ if(i == 5) break print(i) } ## [1] 1 ## [1] 2 ## [1] 3 ## [1] 4
break是退出当前的循环，这里在第5次迭代的时候，判断满足i==5的条件了，于是退出循环，后边的print(i)没有被执行。
for(i in 1:10){ if(i == 5) next print(i) } ## [1] 1 ## [1] 2 ## [1] 3 ## [1] 4 ## [1] 6 ## [1] 7 ## [1] 8 ## [1] 9 ## [1] 10
next是结束当前的迭代，进行下一次迭代。这里在第5次迭代的时候，发现满足了i==5，于是next到下一次迭代，第五次迭代中在next语句之后的print(i)并没有执行。
for(i in 1:3){ for(j in 1:3){ if(i == 2) break print(sprintf("i is %d, j is %d", i, j)) } } ## [1] "i is 1, j is 1" ## [1] "i is 1, j is 2" ## [1] "i is 1, j is 3" ## [1] "i is 3, j is 1" ## [1] "i is 3, j is 2" ## [1] "i is 3, j is 3"
break只影响直接包含它的循环，这里当i==2的时候，从里边的for循环break出来，外边的for依然在继续执行。
for(i in 1:3){ if(i == 2) break for(j in 1:3){ print(sprintf("i is %d, j is %d", i, j)) } } ## [1] "i is 1, j is 1" ## [1] "i is 1, j is 2" ## [1] "i is 1, j is 3"
这里break了外边部的for循环，break之后里边的for循环也不再执行。
for(i in 1:3){ if(i == 2) next for(j in 1:3){ print(sprintf("i is %d, j is %d", i, j)) } } ## [1] "i is 1, j is 1" ## [1] "i is 1, j is 2" ## [1] "i is 1, j is 3" ## [1] "i is 3, j is 1" ## [1] "i is 3, j is 2" ## [1] "i is 3, j is 3"
next在这里只要i==2没有执行，后边i==3的时候得到执行。
R编程
编程就是将一系列的指令组合在一起，告诉电脑做什么，通常会用到上边学到的基本的函数以及流程控制。 我们前边已经学习了R的语法和词语，需要把它们组合，才能够实现丰富多彩的任务。这里展示一些例子。 同样的语法框架下，一个意思有不同的表达方法，使用不同的说法效果不同，这就造就了语言大师和语拙者之间的区别。
例子：仅用加减乘除求X的N次方
仅用加减乘除，求X的N次方。我们先从求10的5次方作为例子来看，我们用加减乘除来求10的五次方，那么我们分5次迭代来实现，如下:
第01次迭代：1*10 第02次迭代：1*10*10 第03次迭代：1*10*10*10 第04次迭代：1*10*10*10*10 第05次迭代：1*10*10*10*10*10
大家可以看出来，这里做了5次迭代，这里边有个规律，就是每次迭代是在上一次迭代的基础上乘以10，总的乘了5次，要知道乘了5次，我们需要一个对象来计数，每次迭代 + 1。当这个迭代次数到达5的时候，我们就停下来，用R语言写出来就是这样子：
# 初始化一个变量，用于随着迭代 + 1，从而实现对迭代次数的记录 die_dai_n <- 0 # 初始化一个变量，用于每次迭代在上一次迭代的基础上 × 10，从而实现5次迭代之后得到10的5次方 n_ci_fang <- 1 # 然后进行重复，到迭代次数die_dai_n 等于 5 的时候停下来 repeat { n_ci_fang <- n_ci_fang * 10 die_dai_n <- die_dai_n + 1 if(die_dai_n == 5) break } # 循环完成之后看看结果 n_ci_fang
R控制台运行效果如下
然后我们推广到，X的N次方，只要把10和5的地方分别换成X和N就可以了，但是要测试先随意定义一个X和N
X <- 2 N <- 3 die_dai_n <- 0 n_ci_fang <- 1 repeat { n_ci_fang <- n_ci_fang * X die_dai_n <- die_dai_n + 1 if(die_dai_n == N) break } n_ci_fang ## [1] 8
我们可以把上边这段代码一般化为一个函数，调用函数修改参数即可方便地求任何X的的任何N次方：
X_N_CiFang <- function(X, N){ die_dai_n <- 0 n_ci_fang <- 1 repeat { n_ci_fang <- n_ci_fang * X die_dai_n <- die_dai_n + 1 if(die_dai_n == N) break } return(n_ci_fang) } X_N_CiFang(X=10, N=5) ## [1] 100000 X_N_CiFang(X=5, N=10) ## [1] 9765625
# 当然，我们也可以用while或者for来实现：
I <- 1; for(i in 1:5) { I <- I*10 } I ## [1] 100000
看看每次迭代都进行了啥：
I <- 1; for(i in 1:5) { I <- I*10 print(sprintf("i = %d, I = %d", i, I)) } ## [1] "i = 1, I = 10" ## [1] "i = 2, I = 100" ## [1] "i = 3, I = 1000" ## [1] "i = 4, I = 10000" ## [1] "i = 5, I = 100000"
如果还不明白，那我们再来让它更清楚一些，引入一个叫迭代(DieDai)的变量，每次迭代这个变量加1：
CiFang <- 1; DieDai <- 0 for(i in 1:5) { DieDai <- DieDai + 1 CiFang <- CiFang * 10 print(sprintf("第%d次迭代，得到CiFang = %d", DieDai, CiFang)) } ## [1] "第1次迭代，得到CiFang = 10" ## [1] "第2次迭代，得到CiFang = 100" ## [1] "第3次迭代，得到CiFang = 1000" ## [1] "第4次迭代，得到CiFang = 10000" ## [1] "第5次迭代，得到CiFang = 100000"
很明显，这里的i是什么并不影响最后的结果，在后边的代码中并不需要用到i的值。
I <- i <- 1 while(i <= 5){ I <- I*10 print(sprintf("i = %d, I = %d", i, I)) i <- i + 1 } ## [1] "i = 1, I = 10" ## [1] "i = 2, I = 100" ## [1] "i = 3, I = 1000" ## [1] "i = 4, I = 10000" ## [1] "i = 5, I = 100000" I <- i <- 1 repeat{ I <- I*10 print(sprintf("i = %d, I = %d", i, I)) if(i == 5) break i <- i + 1 } ## [1] "i = 1, I = 10" ## [1] "i = 2, I = 100" ## [1] "i = 3, I = 1000" ## [1] "i = 4, I = 10000" ## [1] "i = 5, I = 100000"
平时我们在求X的N次方的时候有现成的运算符可以用，表示为X^N即可。如上这样用加减乘除是为了让大家学习循环，锻炼编程思维。
例子：仅用加减乘除求N的阶乘
对于正整数N，有N的阶乘，即
、
、
、
。
比如求5的阶层，我们可以用循环分为五次迭代，如下：
第01次迭代：1*1 第02次迭代：1*1*2 第03次迭代：1*1*2*3 第04次迭代：1*1*2*3*4 第05次迭代：1*1*2*3*4*5
大家可以看出来，这里做了5次迭代，这里边有个规律，就是每次迭代是在上一次迭代的基础上乘以一个数字，这个数字每次迭代是在上一次数字的基础上+1。当这个数字为5的时候就停下来了。 我们用R语言写出来就是这样子。
# 首先定义一个对象，这个对象初始值是0，每次迭代在上一次迭代得到的值的基础上 + 1 die_dai_n <- 0 # 首先定义一个对象，这个对象初始值是1，每次迭代它在上一次迭代的值的基础上 * die_dai_n jie_cheng <- 1 # 然后用repeat进行重复，满足迭代次数 (die_dai_n) 是5的时候停下来，也就是break repeat { die_dai_n <- die_dai_n + 1 jie_cheng <- jie_cheng * die_dai_n if(die_dai_n == 5) break } # 循环结束之后我们看看jie_cheng的值对不对 jie_cheng
看看R控制台运行的效果：
然后我们拓展成求N的阶乘，比如N=10，我们用N替换上边5的位置，把N在前边定义，这样每次可以在前边对N赋值之后运行，即可：
N <- 10 die_dai_n <- 0 jie_cheng <- 1 repeat { die_dai_n <- die_dai_n + 1 jie_cheng <- jie_cheng * die_dai_n if(die_dai_n == N) break } jie_cheng
然后我们可以把这个这个代码放到一个函数里边，方便调用求不同N的阶乘，就是把N当做参数即可：
JieCheng <- function(N){ die_dai_n <- 0 jie_cheng <- 1 repeat { die_dai_n <- die_dai_n + 1 jie_cheng <- jie_cheng * die_dai_n if(die_dai_n == N) break } return(jie_cheng) } JieCheng(5) ## [1] 120 JieCheng(10) ## [1] 3628800
除了用repeat以外，我们可以用while来写，也是一样的，比如求10的阶乘（10!）：
N <- 10 die_dai_n <- 0 jie_cheng <- 1 while(die_dai_n < N) { die_dai_n <- die_dai_n + 1 jie_cheng <- jie_cheng * die_dai_n } jie_cheng ## [1] 3628800
注意，先die_dai_n + 1还是先jie_cheng * die_dai_n的顺序是有影响的，大家可以看看如下的写法，看看循环次数的区别：
N <- 10 die_dai_n <- 1 jie_cheng <- 1 while(die_dai_n <= N) { jie_cheng <- jie_cheng * die_dai_n die_dai_n <- die_dai_n + 1 } jie_cheng ## [1] 3628800
另外还能用for来实现。
jie_cheng <- 1 for(N in 1:10) jie_cheng <- jie_cheng * N jie_cheng ## [1] 3628800
当然，求阶乘还能使用prod或者factorial函数，这里用加减乘除为的是编程示例，现实编程中使用现成的函数可以更高效。一方面，现有的函数经过优化，执行效率更高，一方面，直接使用现成的函数也能节约自己的时间。
例子：通过泰勒公式求e的x次方
根据泰勒公式：
要算出这个问题，我们每次迭代加一项，每一项可以通过其前一项计算得到，每一项在其上项
的基础上，在分子上乘以x，在分母上乘以n，也就是
，而n通过上一迭代的n加+1得到。算到一定的精度就停下来，比如
的绝对值比10*.Machine$double.eps要小即停下来。
edxcf <- function(x=1, delta=10*.Machine$double.eps){ ex <- an <- n <- 1 repeat{ an <- an*(x/n) ex <- ex + an if(abs(an) < delta) break n <- n+1 } return(ex) } options(digits=22) edxcf(1) ## [1] 2.718281828459045534885 edxcf(2) ## [1] 7.389056098930648630585 edxcf(-1) ## [1] 0.3678794411714424450466
例子：斐波拉契数列
比如这里一个求斐波拉契数列(0, 1, 1, 2, 3, 5, 8, 13, 21, …… ,每项是前边两项之和)中不超过某个数(比如1000)的最大的数：
a <- 0; b <- 1 while(a <= 1000){ print(a) c <- a + b a <- b b <- c } ## [1] 0 ## [1] 1 ## [1] 1 ## [1] 2 ## [1] 3 ## [1] 5 ## [1] 8 ## [1] 13 ## [1] 21 ## [1] 34 ## [1] 55 ## [1] 89 ## [1] 144 ## [1] 233 ## [1] 377 ## [1] 610 ## [1] 987
在这段代码的运行中，a是我们要求的，初始值a1是0，初始的时候它的后一位b1是1，这个我们是知道的，于是设定这样的初始值。while后边的第一次运行的时候，我们把a打印到屏幕，也就是打印第一个元素，然后我们求出它的后一位的后一位，这里用c表示，c1=a1+b1。然后往后挪一位，a2=b1，b2=c1，这样就结束了第一个循环，接下来，第二个循环，首先打印到屏幕的a是a2，然后求出c2，然后往后挪一位得到a3，b3。这样到第三次循环的时候，打印出a3，如此往复，直到a大于1000的时候就不再运行了，也就是不打印出a来了。
另外，可以把这些数字放到一个向量F里边
F <- c(0, 1) while(F[length(F)] < 1000){ F <- c(F, F[length(F)] + F[length(F)-1]) } print(F) ## [1] 0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597
取最后一个元素可以使用tail函数，但是似乎要更慢一些。一般轻量运用不用太考虑速度的问题，也可以写成这样：
F <- c(0, 1) while(tail(F, n=1) < 1000){ F <- c(F, sum(tail(F, n=2))) } print(F) ## [1] 0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597
例子：判断一个数是否是质数
如果大于1的自然数N，如果它不能被除了1和它自身以外的其它数整除，也即对于从2到N-1的整数n，都有N除以n的余数不为0，即 N%%n != 0 ，则N是质数，2是最小的质数。而在检验中我们不需要检验那么多，只需要检验2到N的平方根这个区间的n就可以了，因为，如果有一个大于sqrt(N)整数能够整除N，那么其乘数必然是小于sqrt(N)的，也就是在检验sqrt(N)之前的数字的时候已经发现了存在一个数n使得N%%n == 0。
N <- 20 N %% 2:sqrt(N) != 0 ## [1] FALSE TRUE FALSE N <- 17 N %% 2:sqrt(N) != 0 ## [1] TRUE TRUE TRUE for(N in 3:20){ if(all(N %% 2:sqrt(N) != 0)) { print(sprintf("%d是质数", N)) } else { print(sprintf("%d不是是质数", N)) } } ## [1] "3是质数" ## [1] "4不是是质数" ## [1] "5是质数" ## [1] "6不是是质数" ## [1] "7是质数" ## [1] "8不是是质数" ## [1] "9不是是质数" ## [1] "10不是是质数" ## [1] "11是质数" ## [1] "12不是是质数" ## [1] "13是质数" ## [1] "14不是是质数" ## [1] "15不是是质数" ## [1] "16不是是质数" ## [1] "17是质数" ## [1] "18不是是质数" ## [1] "19是质数" ## [1] "20不是是质数"
对于比较大的数来说，我们不需要检验到其平方根那么多，只要找到第一个能够整除它的数就可以停止了。
N <- 20 X <- 1 n <- 2 while(n <= sqrt(N)){ if(N %% n == 0) { X <- 0 break } n <- n + 1 } if(X == 0) { print("N不是质数") } else { print("N是质数") } ## [1] "N不是质数"
while循环运行到找到能够整除N的n就停下来，如果找到了停下来，那么就会把X值改为0，而如果是运行结束了还没有找到能够整除N的n，这时候X的值没有改变，依然是1。 因此，while循环之后通过X的值来判断N是否是质数。
因此，可以写个函数，判断任意一个自然数N是否是质数。
ZS <- function(N, X=1, n=2){ while(n <= sqrt(N)){ if(N %% n == 0) { X <- 0 break } n <- n + 1 } if(X == 0) { print(sprintf("%d非质数", N)) } else { print(sprintf("%d是质数", N)) } } for(N in 10:20) ZS(N) ## [1] "10非质数" ## [1] "11是质数" ## [1] "12非质数" ## [1] "13是质数" ## [1] "14非质数" ## [1] "15非质数" ## [1] "16非质数" ## [1] "17是质数" ## [1] "18非质数" ## [1] "19是质数" ## [1] "20非质数"
例子：取出一个整数各位上的数字
这里实际上用到取整和整除，一个数除以10得到的余数就是最后一位的数字，除以10得到的整数部分就是去掉最后一位数字的数，再一次求余得到十位数的数，以此类推，就能够得到每一位上的数,最后整除得到了0则停止。
N <- 789 N ## [1] 789 N %% 10 ## [1] 9 N <- N %/% 10 N ## [1] 78 N %%10 ## [1] 8 N <- N %/% 10 N ## [1] 7 N %% 10 ## [1] 7 N <- N %/% 10 N ## [1] 0
这里以N为789为例，第一次求余，得到个位上的数；然后整除以10，得到去掉个位上的数剩下的数字，再第二次求余，十位上的数；然后整除以10，去掉十位上的数；然后再除以10求余，得到百位上的数；这时候在整除以10，得到0，说明去掉百位上的数字之后没有数字剩下了，也就停下来了。如此就把一个数字的各位的数字取出来了。
我们把上边的这些代码用while循环来写，仔细想想每一次迭代都做了什么：
N <- 789 while(N > 0){ print(N %% 10) N <- N %/% 10 } ## [1] 9 ## [1] 8 ## [1] 7
我们可以推广，那就是做成一个函数，这样对任意数字都可以方便的进行拆解：
SplitNum <- function(N){ while(N > 0){ print(N %% 10) N <- N %/% 10 } } SplitNum(56789) ## [1] 9 ## [1] 8 ## [1] 7 ## [1] 6 ## [1] 5
我们可以把得到的这些数字放到一个向量里边：
SplitNum <- function(N){ Out <- c() while(N > 0){ Out <- c(N %% 10, Out) N <- N %/% 10 } return(Out) } SplitNum(56789) ## [1] 5 6 7 8 9
先创建一个空的向量Out，每次迭代得到的数字附加到Out的前边 ，最后再把循环完成之后得到的结果作为整个函数运行的结果返回。
如果我要的不是得到每一位上的数字，而是希望对每一位数字进行相加，那么如下：
SumNum <- function(N){ Out <- 0 while(N > 0){ Out <- Out + N %% 10 N <- N %/% 10 } return(Out) } SumNum(56789) ## [1] 35
除了上边的方法以为，可以先将数字转成字符，在是strsplit进行拆分，然后转为数字，然后求和，但是要慢一些，试试如下两行代码可以看出，转为字符拆分的方法要慢两三倍。
system.time(for(i in 1:100000) SumNum(5678)) ## 用户 系统 流逝 ## 0.09 0.00 0.19 system.time(for(i in 1:100000) sum(as.numeric(strsplit(as.character(5678), "")[[1]]))) ## 用户 系统 流逝 ## 0.27 0.00 0.50
除了使用循环以外，还能用向量化的运算：
789 %/% c(100, 10, 1) ## [1] 7 78 789 (789 %/% c(100, 10, 1)) %% 10 ## [1] 7 8 9
向量化整除， 对于3位数，整除的是10^(2:0)，然后再求余，除以十
这是需要知道这个数字的位数的，不知道位数，对任意的正整数，那么可以使用求对数来得到，这样整合起来，做一个函数：
SumNum2 <- function(N){ n <- log(N, 10) %/% 1 + 1 d <- 10^((n-1):0) sum((N %/% d) %/% 10) } SumNum3 <- function(N){ sum((N %/% 10^((log(N, 10) %/% 1):0)) %/% 10) } SumNum2(789) ## [1] 85 SumNum3(789) ## [1] 85
同样可以比较一下这种方法的效率
system.time(for(i in 1:100000000) SumNum(5678)) ## 用户 系统 流逝 ## 111.90 0.19 164.34 system.time(for(i in 1:100000000) SumNum2(5678)) ## 用户 系统 流逝 ## 101.15 0.13 147.14 system.time(for(i in 1:100000000) SumNum3(5678)) ## 用户 系统 流逝 ## 88.14 0.12 127.79
例子：将十进制数转为十六进制数
我们通常写一个数字如果不说明默认是十进制的，也就是456=4×102 + 5×101 + 6×100。而十六进制的45616=4×162 + 5×161 + 6×160。正如我们要求取一个十进制数字的各位数字一样，十进制的数转成16进制的时候，也是一个不断求余和整除的过程。对于16进制数，0~9还不够，后边接着用A-F表示10-15。
为了方便理解，我们先来看看怎么取出十进制数的各位上边的数。
N <- 456 while(N > 0){ print(N %% 10) N <- N %/% 10 } ## [1] 6 ## [1] 5 ## [1] 4
我们来分布展示计算结果，以便大家知道运行过程：
N <- 456 i <- 0 while(N > 0){ i <- i + 1 M <- N %% 10 N <- N %/% 10 print(sprintf("第%d次迭代，除以10取余数得%d，整除以10之后得到%d", i, M, N)) } [1] "第1次迭代，除以10取余数得6，整除以10之后得到45" [1] "第2次迭代，除以10取余数得5，整除以10之后得到4" [1] "第3次迭代，除以10取余数得4，整除以10之后得到0"
同样得道理，在十进制转十六进制时，以16作为除数，这样来求十六进制各位上的数
N <- 2510 while(N > 0){ print(N %% 16) N <- N %/% 16 } ## [1] 14 ## [1] 12 ## [1] 9
这里的2510是十进制数，这里依次得到14、12、9，对应16进制的数字未E、C、9，正好9×16^2+12×16^1+14×16^0 = 2510。
N <- 2510 i <- 0 while(N > 0){ i <- i + 1 M <- N %% 16 N <- N %/% 16 print(sprintf("第%d次迭代，除以16取余数得%d，整除以16之后得到%d", i, M, N)) } ## [1] "第1次迭代，除以16取余数得14，整除以16之后得到156" ## [1] "第2次迭代，除以16取余数得12，整除以16之后得到9" ## [1] "第3次迭代，除以16取余数得9，整除以16之后得到0"
从上边计算的记过，就可以知道，十进制数2510表示为十六进制则为9CE。我们只需要把while循环中得到的余数转成十六进制表示，然后每一次迭代得到的加在前边即可。这样，我们可以写一个函数来转：
Dec2Hex <- function(N){ D2H <- c(0:9, LETTERS[1:6]) names(D2H) <- as.character(0:15) HEX <- "" while(N > 0){ HEX <- paste(D2H[as.character(N %% 16)], HEX, sep="") N <- N %/% 16 } return(HEX) } Dec2Hex(2510) ## [1] "9CE" Dec2Hex(3721) ## [1] "E89" Dec2Hex(255) ## [1] "FF"
同样的原理，我们很容易就能够把十进制的转成其它任何进制的。比如转为2进制则这样：
Dec2Bin <- function(N){ BIN <- "" while(N > 0){ BIN <- paste(N %% 2, BIN, sep="") N <- N %/% 2 } return(BIN) } Dec2Bin(2510) ## [1] "100111001110" Dec2Bin(2) ## [1] "10"
示例：求平方根
比如，我们现在只会算加减乘除，那么怎么计算数字A的平方根呢？ 其中一个方法就是迭代计算，一步一步逼近我们要的值，直到差距小到我们可以接受的范围：
1. 设A的平方根是a
2. 先给定任意一个数X1：如果X1 > a，那么 X1 * X1 > A，而 A/X1 < a， a就介于X1和A/X1之间，X1和A/X1的中点就会更接近a; 同样地，当X1 < a时也是X1和A/X1的中点就会更接近a
3. 于是就有X2 = (X1 + A/X1)/2比X1更接近于a
4. 同样就有X3 = (X2 + A/X2)/2比X2更接近于a
5. 依次不断迭代，计算过程得到数列 Xi+1 = (Xi + A/Xi)/2， 这样, 给定任意一个X1，迭代计算得到X2、X3、X4…，最终会得到无限接近a的值，在达到足够的精度时停下来即可。
下边我们来看看开5的平方，每次都把迭代计算的结果打印出来：
A <- 5 x <- 1 while(TRUE){ x <- (x+A/x)/2 print(x) if(abs(x*x - A) < 0.0000001) { break } } ## [1] 3 ## [1] 2.333333 ## [1] 2.238095 ## [1] 2.236069 ## [1] 2.236068
while(TRUE) 可以换成 repeat，abs()函数求绝对值，当满足条件之后终止循环。
除了上边满足条件的时候break以外，还没有在不满足条件的时候就继续运行：
x <- 1 while(abs(x*x - A) >= 0.0000001){ x <- (x+A/x)/2 print(x) } ## [1] 3 ## [1] 2.333333 ## [1] 2.238095 ## [1] 2.236069 ## [1] 2.236068
很显然，对于任何这样的数列，我们就可以使用迭代法从第一项的值求出各项的值。
例子：水仙花数
如果一个N位数等于它各位上的数字的N次方的和，那么就称这个数就水仙花数。对于N比较小的时候，用R来算就比较方便。如下：
sxhs <- function(N) { for(i in 10^(N-1):(10^N-1)){ n <- 1:N I <- (i %/% 10^(N-n)) %% 10 if(sum(I^N) == i) print(i) } } sxhs(3) ## [1] 153 ## [1] 370 ## [1] 371 ## [1] 407 sxhs(4) ## [1] 1634 ## [1] 8208 ## [1] 9474 sxhs(5) ## [1] 54748 ## [1] 92727 ## [1] 93084
N大于6的时候所需要的时间就太多了，R的for循环还是比较慢的。
sxhs <- function(N) { i <- 10^(N-1) while(i < 10^N){ n <- 1:N I <- (i %/% 10^(N-n)) %% 10 if(sum(I^N) == i) print(i) i <- i + 1 } }
sxhs <- function(N){ i <- 10^(N-1) while(i < 10^N){ j <- 1 S <- 0 I <- i while(j <= N){ S <- S + (I %% 10)^N I <- I %/% 10 j <- j + 1 } if(S == i) print(i) i <- i + 1 } }
都是当N>6就需要花很多时间了。 对于类似的，用R来执行比较慢的，可以考虑使用C++来编写，通过Rcpp来编译和调用。使用?? Rcpp了解更多。
示例：身份证号码校验
根据国家标准(《公民身份证号码》 (GB 11643-1999) ),身份证号从左向右数的第18位为校验码。校验码的计算可以按照如下步骤进行：
1. 前边17位数字ai与其位置的加权因子Wi分别相乘 (从左到右Wi依次为7、9、10、5、8、4、2、1、6、3、7、9、10、5、8、4、2 )
2. 上一步得到的17个数相加，求∑aiWi
3. 上一步得到的数除以11得到余数
4. 用12减去上一步得到的数
5. 上一步得到的数除以11求余数，如果得到10则用X表示
用R来做这件事情，则可以这样：
HM <- "11010519491231002X" W <- c(7, 9, 10, 5, 8, 4, 2, 1, 6, 3, 7, 9, 10, 5, 8, 4, 2) a <- as.numeric(strsplit(as.character(HM), "")[[1]][1:17]) JYM <- (12 - sum(a*W) %% 11) %% 11 if(JYM == 10) JYM <- "X" JYM ## [1] "X"
以上可以做成一个函数，这样方便校验，输入任意身份证好则可以判断最后一位校验码是否复合规则
JY <- function(x){ W <- c(7, 9, 10, 5, 8, 4, 2, 1, 6, 3, 7, 9, 10, 5, 8, 4, 2) A <- strsplit(as.character(x), "")[[1]] if(length(A) != 18) { print("身份证号码位数不对") break } JYM <- as.character((12 - sum(as.numeric(A[1:17])*W) %% 11) %% 11) if(JYM == "10") JYM <- "X" if(JYM == A[18]){ print("身份证号码校验通过") } else { print(sprintf("身份证号码校验码不正确，最后一位应该是:%s。", JYM)) } } JY("11010519491231002X") ## [1] "身份证号码校验通过" JY("11010519491231004X") ## [1] "身份证号码校验码不正确，最后一位应该是:6。"
示例：随机加减算式
小学一年级练习加减法，可能需要一些随机的算式，比如加和不超过某个数的各种数字的加法，随机出两个数相加：
Jia <- function(N){ repeat{ readline() x <- sample(0:N, 1) y <- sample(0:(N-x), 1) cat(sprintf("%d + %d = ", x, y)) }} Jia(N=10)
运行函数之后，每次回车出一个算式。可以把字号调大，这样每次回车屏幕上就显示一个算式。
Jian <- function(N){ repeat{ readline() x <- sample(0:N, 1) y <- sample(0:x, 1) cat(sprintf("%d - %d = ", x, y)) }} Jian(10)
文件与文件夹
如果不指明绝对路径，那么R读和写文件都是在当前的工作目录（Working Directory）下进行的。对于Windows系统，绝对路径就是从磁盘符号开始的，而对于Linux系统，绝对路径则是从根目录/开始的。文件夹及其子文件夹以/分隔，Windows下还能用两个反斜杠\\来分隔。当前工作目录则可以用.代表, 而上一层目录这用..代表，../.. 就是上一级的上一级，以此类推。
创建文件夹与设置工作目录
dir.create函数创建文件夹：
for(i in c("A", "B", "C", "D")){ dir.create(sprintf("D:/dir_not_once_existed/%s", i), recursive = TRUE) }
这里在D盘创建了文件夹dir_not_once_existed，并在其中创建了A、B、C、D四个子文件夹。recursive = TRUE表示可以同时创建多层文件夹。
dir函数查看文件夹下的文件夹及文件：
dir("D:/dir_not_once_existed") ## [1] "A" "B" "C" "D" dir("C:", pattern="Program") ## [1] "Program Files" "Program Files (x86)" "ProgramData"
dir查看指定文件夹中的文件夹或文件，设置pattern可以只查看匹配相应模式的项目，这里看到C盘中名称中包含了Program字符串的。
setwd函数设置工作目录，而getwd查看当前的工作目录：
setwd("D:/dir_not_once_existed/") getwd() ## [1] "D:/dir_not_once_existed" dir() ## [1] "A" "B" "C" "D" setwd("D:/dir_not_once_existed/D") getwd() ## [1] "D:/dir_not_once_existed/D" for(i in c("D1", "D2", "D3", "D4", "D5", "D6")){ dir.create(i) } dir() ## [1] "D1" "D2" "D3" "D4" "D5" "D6"
dir.exists判断文件夹是否存在
dir.exists("D:/dir_not_once_existed/D/D5") ## [1] TRUE dir.exists("D:/dir_not_once_existed/D/D7") ## [1] FALSE
路径指定中，用.表示当前目录，而..表示上一级目录：
setwd(".") getwd() ## [1] "D:/dir_not_once_existed/D" setwd("..") getwd() ## [1] "D:/dir_not_once_existed" setwd("../dir_not_once_existed/D/D5") getwd() ## [1] "D:/dir_not_once_existed/D/D5"
可以用choose.dir()打开对话框选择文件夹，比如setwd(choose.dir())通过对话框选择文件夹的形式设置工作目录。同样的，有file.choose()用于选择单个文件，而choose.files()选取多个文件。
简单文件操作
file.create创建空的文件，cat往文件里写入内容，file.show查看文件内容，file.info查看文件的信息，file.remove删除文件，file.rename重命名文件，file.copy复制文件，file.append将后边文件的内容附加到前边的文件后边。
setwd("D:/dir_not_once_existed/D") dir() ## [1] "D1" "D2" "D3" "D4" "D5" "D6" file.create(c("A", "B", "C")) ## [1] TRUE TRUE TRUE dir() ## [1] "A" "B" "C" "D1" "D2" "D3" "D4" "D5" "D6" file.info(dir()) ## size isdir mode mtime ctime atime exe ## A 0 FALSE 666 2022-12-11 20:59:35 2022-12-11 20:59:35 2022-12-11 20:59:35 no ## B 0 FALSE 666 2022-12-11 20:59:35 2022-12-11 20:59:35 2022-12-11 20:59:35 no ## C 0 FALSE 666 2022-12-11 20:59:36 2022-12-11 20:59:36 2022-12-11 20:59:36 no ## D1 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no ## D2 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no ## D3 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no ## D4 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no ## D5 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no ## D6 0 TRUE 777 2022-12-11 20:49:02 2022-12-11 20:49:02 2022-12-11 20:49:02 no file.remove(c("C")) ## [1] TRUE dir() ## [1] "A" "B" "D1" "D2" "D3" "D4" "D5" "D6" file.rename(c("A", "B"), c("X", "Y")) ## [1] TRUE TRUE dir() ## [1] "D1" "D2" "D3" "D4" "D5" "D6" "X" "Y" cat("X\tY\n", file = "A") cat("BBBB\t4\n", file = "B") cat("C\t10\n", file = "C") cat("DD\t5\n", file = "D") file.show(c("A", "B")) ## 这里会有两个弹窗，展示文件A和B的内容 file.append("A", "B") file.show(c("A", "B")) ## 这里会有两个弹窗，展示文件A和B的内容，留意A后的后边附加了B的内容 file.append("A", c("C", "D")) file.show("A") ## 这里会有一个弹窗，展示文件A的内容，留意A后的后边依次附加了C和D的内容 file.copy("A", "../F") ## [1] TRUE dir("D:/dir_not_once_existed/") ## [1] "A" "B" "C" "D" "F"
文件重命名, 将A重命名为X，B重命名为Y，一一对应的。在interactive模式下用R的时候，可以用file.show这个命令来展示文件的内容。
文件读写之readLines
readLines函数用于逐行读入纯文本的文件，得到一个向量，每行是一个元素
setwd("D:/dir_not_once_existed/D") X <- readLines("A") X ## [1] "X\tY" "BBBB\t4" "C\t10" "DD\t5" str(X) ## chr [1:4] "X\tY" "BBBB\t4" "C\t10" "DD\t5"
文件读写之read.table
read.table序列函数用于读入表格，设置sep参数指定列的分隔符。read.delim默认使用制表符分隔，在读入制表符分隔的表格的时候，使用read.delim更方便。
Y <- read.delim("A") Y ## X Y ## 1 BBBB 4 ## 2 C 10 ## 3 DD 5 str(Y) ## 'data.frame': 3 obs. of 2 variables: ## $ X: chr "BBBB" "C" "DD" ## $ Y: int 4 10 5 Z <- read.delim("A", sep=";", header=FALSE, stringsAsFactors=TRUE) Z ## V1 ## 1 X\tY ## 2 BBBB\t4 ## 3 C\t10 ## 4 DD\t5 str(Z) ## 'data.frame': 4 obs. of 1 variable: ## $ V1: Factor w/ 4 levels "BBBB\t4","C\t10",..: 4 1 2 3
read.delim 默认情况下第一行当作表头，当没有表头的时候需要设置header=FALSE。stringsAsFactors一般设置成为FALSE，这个可以在一个大脚本之初就设置好，默认与全局设置一致。 read.delim默认是以制表符（\t）分隔的，设成了;之后整个文件就只有1列了，可以看到制表符读入为\t。
read.delim属于read.table系列的，还有read.csv、 read.delim2、read.csv2等，各个函数的默认参数不同，对不同的文件使用起来的方便程度不同。 read.csv适合于与csv文件，即以逗号分隔的文件。
文件读写之write.table
write.table用于写出表格。
Z <- data.frame(X=letters[1:5], Y=1:5) Z ## X Y ## 1 a 1 ## 2 b 2 ## 3 c 3 ## 4 d 4 ## 5 e 5 write.table(Z, "Z.tsv", sep="\t", quote=FALSE, row.names=FALSE) file.show("Z.tsv") ## 这里将会弹出一个框展示Z的内容
write.table有几个参数需要知道：一个是sep=“\t”指使用制表符分隔，有的系统默认是空格，需要设定清楚，可以根据需要使用别的，比如逗号等等的；quote=FALSE是指字符写出的时候是不加引号，一般是不加引号的，需要加则改为TRUE；row.names=FALSE是指不把行名写出，如果写出，通常会导致第一行向前错位一列。 可以尝试不同的参数看看效果。
在Windows系统中，write.table写出的换行符是“\r\n”，有时候写出的文件需要放到Linux系统中使用，这时候换行符就需要用“\n”，这时候简单设置eol=“\n”是不能够实现目的的，需要以二进制的方式打开文件写入,这样：
File <- file("Zn.tsv", "wb") write.table(Z, File, sep="\t", quote=FALSE, row.names=FALSE) close(File)
用Notepad++打开文件Zn.tsv在下边信息栏的右边就能看到Unix(LF)的字样了，而上边的Z.tsv文件则显示为Windows(CR LF)。CR是回车(\r)，LF是换行(\n)，Windows采用CR加LF表示换行，而Linux使用LF表示换行，而MAC机则是使用CR表示换行。 wb表示以二进制写方式打开，只能写， 文件不存在则创建，文件存在则覆盖。
write.table之后，要close，不然文件大小显示为0，也打不开。
文件读写之writeLines
writeLines将一个向量的元素逐个写入文件，每个元素一行。
writeLines(Z$X, "Z.txt") file.show("Z.txt")
关于文件换行的问题，参照write.table的做法。
超大文件读写之data.table::fread与data.table::fwrite
# 这里我们先写出一个测试文件，1亿行，2列 set.seed(2025); Text <- sample(letters, 10^8, replace=TRUE) set.seed(5202); Num <- sample(1:10^9, 10^8) dat <- data.frame(Text, Num) # 把它写到硬盘，测试的时候注意写入的目录所在硬盘要有足够的空间，注意这个文件名，是否和自己已有的文件名有重名的，以免覆盖重要文件 fwrite(dat, "test_big_table.tsv", row.names=FALSE, sep="\t") # 可以看到这个文件有 file.info("test_big_table.tsv") ## size isdir mode mtime ctime atime exe uname udomain ## test_big_table.tsv 1288885892 FALSE 666 2025-10-31 12:02:06 2025-10-31 10:05:20 2025-10-31 no Everyone Dat <- fread("test_big_table.tsv") # 可以看到，读入的数据除了data.frame这个class以外，还有data.table这个class str(Dat)
自己查看说明文档去了解更多的参数细节。
关于速度的测试看这里，同学们可以不用去亲自测试，了解一下这个速度差别即可。不同的电脑所消耗的时间是不一样的。由于文件比较大，我在测试的过程中把原先内存里的dat先删除了，并且gc()做了内存整理再读入，以免超内存报错。详见如下记录：
library(data.table) setwd("E:/test") set.seed(2025); Text <- sample(letters, 10^8, replace=TRUE) set.seed(5202); Num <- sample(1:10^9, 10^8) dat <- data.frame(Text, Num) str(dat) ## 'data.frame': 100000000 obs. of 2 variables: ## $ Text: chr "m" "l" "d" "z" ... ## $ Num : int 403315539 862071504 672821329 482661279 835304693 301825546 731876629 64240564 46055120 788547792 ... system.time(write.table(dat, "test_01.tsv", row.names=FALSE, sep="\t", quote=FALSE)) ## 用户 系统 流逝 ## 315.47 56.78 509.04 system.time(fwrite(dat, "test_02.tsv", row.names=FALSE, sep="\t")) ## ## Written 61.3% of 100000000 rows in 14 secs using 6 threads. maxBuffUsed=40%. ETA 9 secs. ## ## ## 用户 系统 流逝 ## 5.40 0.78 22.41 file.info(c("test_01.tsv", "test_02.tsv")) ## size isdir mode mtime ctime atime exe ## test_01.tsv 1288885892 FALSE 666 2025-10-31 12:02:06 2025-10-31 10:05:20 2025-10-31 no ## test_02.tsv 1288885892 FALSE 666 2025-10-31 12:02:22 2025-10-31 10:12:43 2025-10-31 no ## uname udomain ## test_01.tsv Everyone ## test_02.tsv Everyone rm(list=ls()) gc() ## used (Mb) gc trigger (Mb) max used (Mb) ## Ncells 369657 19.8 70937852 3788.5 100370394 5360.4 ## Vcells 67790997 517.3 725147016 5532.5 728932689 5561.4 system.time(dat01 <- read.delim("test_01.tsv")) ## 用户 系统 流逝 ## 286.95 8.55 300.43 system.time(dat02 <- fread("test_01.tsv")) ## |--------------------------------------------------| ## |==================================================| ## 用户 系统 流逝 ## 6.69 2.17 38.33 str(dat01) ## 'data.frame': 100000000 obs. of 2 variables: ## $ Text: chr "m" "l" "d" "z" ... ## $ Num : int 403315539 862071504 672821329 482661279 835304693 301825546 731876629 64240564 46055120 788547792 ... str(dat02) ## Classes ‘data.table’ and 'data.frame': 100000000 obs. of 2 variables: ## $ Text: chr "m" "l" "d" "z" ... ## $ Num : int 403315539 862071504 672821329 482661279 835304693 301825546 731876629 64240564 46055120 788547792 ... ## - attr(*, ".internal.selfref")=<externalptr>
文件读写之图片保存
pdf打开设备作图保存为pdf文件，结束之后，用dev.off()关闭设备。
pdf("plot.pdf", height=7, width=7, family="GB1") plot(x=1:10, y=10:1, xlab="横坐标标题", ylab="纵坐标标题") dev.off()
对于图片中含有中文字符的时候，这里需要加上family=“GB1”这个参数，不然中文无法显示。 除了pdf以外，还有png，jepg，tif等等图片格式，也是同样的方法。从pdf()那一行开始，到dev.off()，整个保存图片的工作才算结束了。相当于pdf打开了一个设备，往后所有的图片都在这个pdf里边画，到这里关闭了这个设备，文件得到最终保存。要注意：最后的dev.off()不能少，不执行这条命令文件大小显示为0，无法打开。
再比如，对于一个ggplot2做的图, 保存图片除了用ggsave以外，还能用pdf这样打开一个设备的方式保存。
library(ggplot2) g <- ggplot() + geom_point(aes(x=1:5, y=1:5)) png("ggplot_save_test_01.png") g dev.off() ggsave("ggplot_save_test_02.png", plot=g)
图片保存在了当前的工作目录，可以使用getwd()查看目前的工作目录。比如这里上边将得到两个png文件。 可以用? png 和 ? ggsave查看更多的可选参数。注意，png和pdf默认的长宽的单位不一样，png默认的单位是像素。
文件读写之save与load
用save保存对象，然后使用load加载：
setwd("D:/dir_not_once_existed/") rm(list=ls()) x <- 1:10 y <- data.frame(A=letters[1:10], B=1:10) ls() ## [1] "x" "y" save(x, y, file="r_objects.RData") dir() ## [1] "A" "B" "C" "D" "F" ## [6] "r_objects.RData" rm(list=ls()) ls() ## character(0) load("r_objects.RData") ls() ## [1] "x" "y"
文件中R代码的调用
这里将cat('大家好\n')这个代码写入文件code_test.r：
setwd("D:/dir_not_once_existed/") cat("cat('大家好\n')", file = "code_test.r") dir() ## [1] "A" "B" "C" "code_test.r" "D" ## [6] "F" "r_objects.RData"
是使用source调用保存在文件中的R代码：
source("code_test.r") ## 大家好
R脚本
除了如上在R控制台运行R的代码以外，还能把代码放到一个文件里，作为脚本的形式运行。
脚本及调用
比如这个traditional_greeting.R脚本：
- traditional_greeting.R
cat("Hello World")
比如把这个脚本写到D盘，其路径为D:/tranditional_greeting.R：
writeLines('cat("Hello World")', "D:/tranditional_greeting.R")
在windows下，使用 Windows + r 调出运行窗口，然后输入cmd，点击确定就可以打开，就可以打开控终端，在终端输入如下命令，即可运行，得到脚本运行的结果，比如这里是把Hello World输出到屏幕上：
D:/Program/R/R-4.0.3/bin/Rscript.exe D:/tranditional_greeting.R
这里的D:/Program/R/R-4.0.3/bin/Rscript.exe是我电脑上的安装路径，相应的要换成大家自己电脑上的路径。如果路径中有空格则将整个路径用引号引起来。这里用绝对路径，这样不管当前工作目录在哪里都可以找到它。
比如这里的运行情况如下：
如果这个路径添加到了PATH环境变量中，则可以直接使用Rscript.exe(.exe可省略)进行,在windows终端CMD中输入如下命令并执行：
setx PATH "%PATH%;D:/Program/R/R-4.0.3/bin/"
重新打开的终端CMD里就可以简单的用Rscript运行脚本了，如果将目录切换到脚本所在的目录，则直接用用脚本文件名即可：
脚本变量传入
另外，还能将一些变量作为参数的形式传入，只需加上args ← commandArgs(), 然后使用args[i]来取出第i个参数。比如试试这个test.R脚本，存在D盘：
- test.R
args <- commandArgs(trailingOnly = FALSE) #trailingOnly默认的就是FALSE for(i in 1:length(args)){ print(sprintf("第%d个参数是：%s", i, args[i])) }
设置trailingOnly = TRUE那么只有–args后边的参数被返回，通常设为TRUE，比较方便使用，看看效果：
- test.R
args <- commandArgs(trailingOnly = TRUE) # for(i in 1:length(args)){ print(sprintf("第%d个参数是：%s", i, args[i])) }
比如上边的开平方的函数，我们可以装入一个脚本，方便调用。
- sqrt.R
args <- commandArgs(trailingOnly = TRUE) A <- as.numeric(args[1]) i <- A/2 while(abs(i*i - A) >= 0.0000001){ # 当没有满足条件的时候，继续，这样把每次迭代的结果都输出 i <- (i+A/i)/2 } cat(sprintf("One sqrt root of %s is %s", A, i)) # 整理一下输出结果
这里将此文件保存在上文中已经设置的工作目录中，则可用于开平方啦！
Rscript sqrt.R 99 Rscript sqrt.R 33
R脚本的编码问题
当R脚本中含有汉字的时候，再Windows中，R脚本文件的编码方式需要是UTF-8、GB2312等支持中文字符的编码方式，这个使用Notepad++编辑脚本文件的时候可以自己设置：菜单栏 → 编码 → 编码字符集 → 中文 → GB2312(简体中文)。一些同学发现脚本无法运行，一些同学可能发现运行之后文件夹名称等等的是乱码，其问题可能就出在这个R脚本的编码问题上，试试改一下编码方式看看问题是否解决。
R内调用R脚本
为了调试的方便，我们可以从R控制台运行R的脚本，这时候用到sys这个包的exec_wait函数，如果没有这个包则先安装一下，如果还没有上文提到的测试脚本，这个里加一下：
if(!requireNamespace("sys")) install.packages("sys") if(!file.exists("tranditional_greeting.R")) writeLines('cat("Hello World")', "tranditional_greeting.R")
假如你已经像上文中那样将Rscript.exe文件所在目录加到了PATH变量中，则可以在R中运行这样的命令：
sys::exec_wait("Rscript tranditional_greeting.R")
能够看到这里把运行的结果输入到了屏幕，最后把程序退出状态(exit code)也输出到了屏幕上，0表示程序运行成功。更多信息可以使用? exec_wait查看更多。
R启动设置
R启动时默认工作路径是我的文档，在我的文档这个文件夹里边添加一个名为.Rprofile的文件，在这个文件里可以设置想要R启动时默认设置的内容。比如改变工作路径。
- .Rprofile
# 设置默认的镜像，比如CRAN的和Bioconductor的，设成国内的速度估计都是足够的，试试，有必要再换成别的 options(repos="https://mirrors.hust.edu.cn/CRAN") options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") # 采用科学计数法显示的位数，默认是数字达到8位的时候就会显示为科学计数 options(scipen=99) # 最大有效数字位数，显示精度，不影响计算精度，可设 0-22 options(digits=22) # 更改默认的工作目录，这里可以更改为你自己设定的，下一行更改后把前边的#去掉才生效，下同 # setwd("D:/") # 设置R包的路径，这个一般不用设置，但是如果系统默认的路径不合适时可以加上 # .libPaths("D:/Program/R/R-4.2.0/library") # 可以默认启动加载一些包，比如经常要用到的包，这样方便， # 但是这样要注意，比如分享代码的时候别忘记把这些包加上，所以我一般都是用到再加载，以便记录在代码记录里 # library(ggplot2) # library(plyr) # library(reshape) # library(pheatmap)
数学运算（常用函数）
先创建四个向量，用于展示示例：
set.seed(2022); (a <- rnorm(10, mean=1)) ## [1] 1.37697212 1.30154837 -0.09802317 -0.13040590 -1.79653432 1.72057350 1.93912102 0.77062225 2.75913135 1.11736679 set.seed(2022); (x <- sample(10)) ## [1] 4 3 6 7 9 10 2 8 5 1 set.seed(2022); (y <- sample(1:100,10)) ## [1] 100 51 78 55 75 68 6 63 98 14 set.seed(2202); (z <- sample(1:100,10)) ## [1] 99 8 74 87 39 45 55 16 96 22
round(x, digits = 0)进行四舍五入到小数点后digits位，要注意的是对5的舍和入可能在不同的系统上不同，比如有的时候用R算出来和用EXCEL算出来就不一样，需要留意，使用合适的函数；signif(x, digits = 6)四舍五入到digits位有效数字；floor(x)得到不比x大的最大整数；ceiling(x)得到不比x小的最小整数；trunc(x)得到像0一侧的离x最近的整数。
round(a, digits = 6) ## [1] 1.376972 1.301548 -0.098023 -0.130406 -1.796534 1.720573 1.939121 0.770622 2.759131 1.117367 signif(a, digits = 6) ## [1] 1.3769700 1.3015500 -0.0980232 -0.1304060 -1.7965300 1.7205700 1.9391200 0.7706220 2.7591300 1.1173700 round(a*10^(0:9), digits = 6) ## [1] 1.376972 13.015484 -9.802317 -130.405904 -17965.343199 172057.349841 1939121.023009 7706222.532925 275913134.696347 1117366786.802848 signif(a*10^(0:9), digits = 6) ## [1] 1.37697 13.01550 -9.80232 -130.40600 -17965.30000 172057.00000 1939120.00000 7706220.00000 275913000.00000 1117370000.00000 round(a) ## [1] 1 1 0 0 -2 2 2 1 3 1 floor(a) ## [1] 1 1 -1 -1 -2 1 1 0 2 1 ceiling(a) ## [1] 2 2 0 0 -1 2 2 1 3 2 trunc(a) ## [1] 1 1 0 0 -1 1 1 0 2 1 round(c(1234, 3456), -2) ## [1] 1200 3500
sum(x)、prod(x)、mean(x)、var(x)、sd(x)、median(x)、max(x)、min(x)用于求数值型向量x的所有元素的总和、乘积、平均值、方差、标准差、中位数、最大值、最小值
sum(x) ## [1] 55 prod(x) ## [1] 3628800 mean(x) ## [1] 5.5 var(x) ## [1] 9.166667 sd(x) ## [1] 3.02765 median(x) ## [1] 5.5 max(x) ## [1] 10 min(x) ## [1] 1
range(x)返回x的最大值和最小值； diff(x)求向量x的元素之差；quantile(x)求x的分位数，可以通过probs参数设置相应的分位，默认是c(0, 0.25, 0.5, 0.75, 1)；rank(x)求x的元素的排序序号；scale(x)将x的值进行标准化，如x是一个矩阵，则是按列进行标准化，得到z-score。
range(x) ## [1] 1 10 diff(x) ## [1] -1 3 1 2 1 -8 6 -3 -4 quantile(x, probs=seq(0, 1, 0.1)) ## 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% ## 1.0 1.9 2.8 3.7 4.6 5.5 6.4 7.3 8.2 9.1 10.0
var(x, y) 或者cov(x, y)、cor(x, y) 用于求x和y之间的协方差、相关性系数。cor默认是求皮尔逊相关系数，可以通过method参数改为求spearman或者kendall系数。
var(x, y) ## [1] 48.55556 cov(x, y) ## [1] 48.55556 cor(x, y) ## [1] 0.5129468 cor(x, y, method = "spearman") ## [1] 0.4181818
var、cov、cor还能用于数据框或者矩阵，这样计算各列两两之间的相关性系数：
set.seed(2018); X <- matrix(sample(1:100), 20, 5, dimnames=list(letters[1:20], LETTERS[1:5])); X ## A B C D E ## a 15 35 67 100 28 ## b 55 43 6 65 7 ## c 3 82 14 46 63 ## d 12 49 91 96 90 ## e 18 42 56 27 62 ## f 85 86 95 8 83 ## g 73 21 57 47 44 ## h 72 61 69 31 71 ## i 32 25 9 98 10 ## j 92 74 87 52 66 ## k 88 84 26 64 4 ## l 79 2 41 80 40 ## m 97 77 34 89 45 ## n 19 99 81 20 5 ## o 22 54 33 59 93 ## p 1 23 39 70 50 ## q 94 53 11 30 36 ## r 78 29 16 75 68 ## s 13 48 17 38 60 ## t 51 76 58 24 37 cor(X) ## A B C D E ## A 1.00000000 0.119992673 0.01330526 -0.1110338 -0.091366317 ## B 0.11999267 1.000000000 0.30348296 -0.5093358 -0.001618744 ## C 0.01330526 0.303482964 1.00000000 -0.2518418 0.344120441 ## D -0.11103385 -0.509335763 -0.25184178 1.0000000 -0.143649964 ## E -0.09136632 -0.001618744 0.34412044 -0.1436500 1.000000000
pmax、pmin是向量化的运算，对给的的两个或者两个以上的向量，依次对应元素秋最大值。
max(x, y, z) ## [1] 100 pmax(x, y, z) ## [1] 100 51 78 87 75 68 55 63 98 22 pmin(x, y, z) ## [1] 4 3 6 7 9 10 2 8 5 1
cumsum和cumprod用于逐元素累加或者累乘。
sum(1:10) ## [1] 55 cumsum(1:10) ## [1] 1 3 6 10 15 21 28 36 45 55 prod(1:10) ## [1] 3628800 cumprod(1:10) ## [1] 1 2 6 24 120 720 5040 40320 362880 3628800
log(x)、log10(x)、log2(x)、log(x, a)分别求x的以e（自然常数）、10、2、a为底的对数，exp(x)求自然常数e的x次方，这些都是向量化的。
log(1:10) ## [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101 2.0794415 2.1972246 2.3025851 log10(1:10) ## [1] 0.0000000 0.3010300 0.4771213 0.6020600 0.6989700 0.7781513 0.8450980 0.9030900 0.9542425 1.0000000 log(1:10, 100) ## [1] 0.0000000 0.1505150 0.2385606 0.3010300 0.3494850 0.3890756 0.4225490 0.4515450 0.4771213 0.5000000
sin(x)、cos(x)、tan(x)、asin(x)、acos(x)、atan(x)求x的正弦、余弦、正切、反正弦、反余弦、反正切。另外还有atan2(y, x)求x轴与原点到点(x, y)的向量之间的夹角，当y和x都是正数的时候，atan2(y, x) == atan(y/x)。 同样，这些函数也是向量化的。
sin(1:10) ## [1] 0.8414710 0.9092974 0.1411200 -0.7568025 -0.9589243 -0.2794155 0.6569866 0.9893582 0.4121185 -0.5440211 asin(seq(-1, 1, length=10)) ## [1] -1.5707963 -0.8911225 -0.5890310 -0.3398369 -0.1113410 0.1113410 0.3398369 0.5890310 0.8911225 1.5707963 atan2(1, 1) ## [1] 0.7853982 atan2(-1, -1) ## [1] -2.356194
字符操作（常用函数）
(x <- "ATTATGCATCGGC") ## [1] "ATTATGCATCGGC" (y <- c("CGATTACGGC", "ATGAATGGC", "ATGTCGC", "ATNAGGGC")) ## [1] "CGATTACGGC" "ATGAATGGC" "ATGTCGC" "ATNAGGGC"
字符计数 (Count the Number of Characters)
nchar函数返回字符串中的字符数，这也是向量化的函数。注意并不是用length，length是得到向量y的长度，即元素个数。
nchar(y) ## [1] 10 9 7 8 length(y) ## [1] 4
提取或改变字符串的部分内容(Substrings of a Character Vector)
使用substr可以提取字符串的部分内容，加上赋值符号可以改变这部分内容。 substr是向量化的函数，对每个元素都进行同样操作，返回对应顺序结果。 substring默认设置了last参数，所以可以不提供，一般来说，字符串没有这么长，所以当从一个位置取到最后一位的时候，用substring更便捷。
substr(y, 2, 3) ## [1] "GA" "TG" "TG" "TN" substr(y, 2:5, 3:4) ## [1] "GA" "GA" "" "" substring(y, 4:7) ## [1] "TTACGGC" "AATGGC" "TCGC" "AGGGC" substring("ABCDEFGH", 4:7) ## [1] "DEFGH" "EFGH" "FGH" "GH" Y <- y substr(Y, 2, 3) <- c("11", "22", "3", "456789") Y ## [1] "C11TTACGGC" "A22AATGGC" "A3GTCGC" "A45AGGGC" Y <- y substring(Y, 6) <- c("11", "22", "33", "456789") Y ## [1] "CGATT11GGC" "ATGAA22GC" "ATGTC33" "ATNAG456"
可以看到，在替换的时候，最多只能替换为同样长度的子字符串。而使用substring进行的时候，被替换的位数最多是替换字符的位数。
字符串分割(Split the Elements of a Character Vector)
strsplit用来将字符串按照指定的分隔字符(串)进行分割，返回结果为列表的形式。
strsplit(x, "") ## [[1]] ## [1] "A" "T" "T" "A" "T" "G" "C" "A" "T" "C" "G" "G" "C" strsplit(x, "T") ## [[1]] ## [1] "A" "" "A" "GCA" "CGGC" strsplit(y, "T") ## [[1]] ## [1] "CGA" "" "ACGGC" ## ## [[2]] ## [1] "A" "GAA" "GGC" ## ## [[3]] ## [1] "A" "G" "CGC" ## ## [[4]] ## [1] "A" "NAGGGC"
模式匹配与替换(Pattern Matching and Replacement)
grep(pattern, x)用于对x进行搜索，看其元素中是否包含正则表达式pattern。返回含pattern的元素的序号(index)。如果加上value=TRUE，返回含pattern的元素。而grepl则是返回逻辑型向量。
grep("CG", y) ## [1] 1 3 grep("CG", y, value = TRUE) ## [1] "CGATTACGGC" "ATGTCGC" grepl("CG", y) ## [1] TRUE FALSE TRUE FALSE grepl("^CG", y) ## [1] TRUE FALSE FALSE FALSE grepl("C.G", y) ## [1] TRUE FALSE FALSE FALSE gsub("CG", "_CG_" , y) ## [1] "_CG_ATTA_CG_GC" "ATGAATGGC" "ATGT_CG_C" "ATNAGGGC"
当整个向量中固定某个模式都有且只有唯一匹配的时候，可以使用gsub提取，但是这时候没有匹配的项目将会被整个返回来，只有非常确定结果的时候才能用，否则还是使用上文中提到的stringr的函数要好一些。比如:
gsub(".*([0-9]{12}).*", "\\1", "小明的学号是202201010001呀") ## [1] "202201010001" gsub(".*([0-9]{12}).*", "\\1", c("小明的学号是2022010100013呀", "是的")) ## [1] "022010100013" "是的"
使用如下命令查看更多的关于正则表达式的内容：
? regex
字符串串联 (Concatenate Strings)
(z <- letters[1:5]) ## [1] "a" "b" "c" "d" "e" paste(x, z, sep="_") #向量化的函数，将两个向量的元素粘在一起，以特定字符相隔 ## [1] "ATTATGCATCGGC_a" "ATTATGCATCGGC_b" "ATTATGCATCGGC_c" "ATTATGCATCGGC_d" "ATTATGCATCGGC_e" paste(z, collapse="_") #设置collapse参数，将向量中所有的元素粘在一起，以特定字符相隔 ## [1] "a_b_c_d_e"
大小写转换
sprintf
比如，现在有一个表格，记录了几个人的成绩，如下：
姓名 数学 语文 英语 张三 99 85 81 李四 90 86 82 王五 91 87 83 赵六 92 98 84
那如果现在我要计算机一个个把这些成绩念出来，说“张三数学99、语文85、英语81”这样，这里边有姓名和三个科目的成绩是可变的，这时候用sprintf来产生这样的格式化的句子：
XM <- c("张三", "李四", "王五", "赵六") SX <- c(99, 90, 91, 92) YW <- c(85, 86, 87, 98) YY <- c(81, 82, 83, 84) sprintf("%s数学%d、语文%d、英语%d", XM, SX, YW, YY) ## [1] "张三数学99、语文85、英语81" "李四数学90、语文86、英语82" "王五数学91、语文87、英语83" "赵六数学92、语文98、英语84"
sprintf函数用于将一些变量格式化并加入字符串中。按照格式相对应位置的变量格式化之后加入字符串。比如：
sprintf("π is %f", pi) ## [1] "π is 3.141593"
将pi的值按照%f规定的格式放到它所在的位置。
sprintf("本课程老师姓%s名%s", "王", "弘") ## [1] "本课程老师姓王名弘"
这里两个%s占据的位置由后边的两个对应顺序的值来替代。
格式以%开头，常用的格式有：
%f代表带小数点的数字，默认小数点后为6位；%.nf代表小数点后保留n位有效数字，比如%.2f代表小数点后保留2位有效数字；
%d代表整数，在相应的位置插入整数；%nd代表整数的最小长度，如果整数不够，默认在前边加空格，比如%2d，就会在0-9的整数前边加上一个空格，而10以上的数字则不在前边加空格；可以在n的前边加上0则前边以0补足数字长度。比如sprintf(”%03d”, 2)就得到“002”
%s代表字符串。
更细致的格式定义可以通过? sprintf查看使用文档了解。
sprintf("%s is %f feet tall\n", "Seven", 7.1) ## [1] "Seven is 7.100000 feet tall\n"
以上代码中，“Seven”替换了%s，7.1替换了%f，替换的时候根据格式7.1后边加了五个0，符合%f默认小数点后6位的格式。
sprintf("%s is %f and %s is %d", "A", 99, "B", 50) ## "A is 99.000000 and B is 50"
如上代码中，%s、%f、%s、%d四个个地方又后边的“A”、99、“B”、50依次按格式替换。
Mean <- 10; SD <- 5 sprintf("%.2f±%.2f(Mean±SD)", Mean, SD) ## [1] "10.00±5.00(Mean±SD)"
以上代码中，Mean的值替代了前边的%.2f，SD的值，替换了后边一个%.2f，都是小数点后保留两位。这和引号里句子中括号中的Mean和SD没有关系，那是写定了的，不会被替换。 这里的变量可以用别的变量名，效果也是一样的，比如：
X <- 10; Y <- 5 sprintf("%.2f±%.2f(Mean±SD)", X, Y) ## [1] "10.00±5.00(Mean±SD)"
A <- 1; B <- 2; C <- 3 sprintf("A is %.2f, B is %02d, C is %s", A, B, C) ## [1] "A is 1.00, B is 02, C is 3"
以下展示不同形式的pi,看看不同的设置的格式化效果：
sprintf("%f", pi) ## [1] "3.141593" sprintf("%.3f", pi) ## [1] "3.142" sprintf("%1.0f", pi) ## [1] "3" sprintf("%5.1f", pi) ## [1] " 3.1" sprintf("%05.1f", pi) ## [1] "003.1" sprintf("%+f", pi) ## [1] "+3.141593" sprintf("% f", pi) ## [1] " 3.141593" sprintf("%-10f", pi) ## [1] "3.141593 " sprintf("%e", pi) ## [1] "3.141593e+00" sprintf("%E", pi) ## [1] "3.141593E+00" sprintf("%g", pi) ## [1] "3.14159" sprintf("%g", 1e6 * pi) ## [1] "3.14159e+06" sprintf("%.9g", 1e6 * pi) ## [1] "3141592.65" sprintf("%G", 1e-6 * pi) ## [1] "3.14159E-06"
sprintf是个向量化的函数，格式对应的向量如果有多个元素，那么就会得到一个向量。如下：
A <- sprintf("%03d", seq(0,100,5)) print(A) ## [1] "000" "005" "010" "015" "020" "025" "030" "035" "040" "045" "050" ## [12] "055" "060" "065" "070" "075" "080" "085" "090" "095" "100" sprintf("x is %d and y is %02d", 1:10, 1:3) ## Error in sprintf("x is %d and y is %02d", x, y) : 不是所有的参数都一样长 sprintf("x is %d and y is %02d", 1:10, 1:5) ## [1] "x is 1 and y is 01" "x is 2 and y is 02" "x is 3 and y is 03" "x is 4 and y is 04" "x is 5 and y is 05" ## [6] "x is 6 and y is 01" "x is 7 and y is 02" "x is 8 and y is 03" "x is 9 and y is 04" "x is 10 and y is 05"
如果要得到%本身，则连用两个，也就是%%，比如：
sprintf("%.2f is %.2f%%", 0.1, 0.1*100) ## [1] "0.10 is 10.00%"
apply序列函数
R里尽量避免for循环，运行速度比较慢。需要对一系列值做同样的操作的时候，可以使用apply序列的函数。
apply
apply其用法如下：
apply(X, MARGIN, FUN, ..., simplify = TRUE)
第一个参数X表示要操作的对象X，MARGIN表示对行（1）或者对列（2）执行后边函数FUN，FUN可以给一个函数名，也可以给一个函数，…的地方可以提供更多的FUN的参数，而simplify=TRUE表示结果在可以简化的时候进行简化，默认是TRUE的，这时候结果得到的是一个向量，而如果改成FALSE之后得到的结果是一个列表。
比如，这里我们构建一个数据库用于演示, 做一个示例数据框，5列，10行。
dat <- data.frame(matrix(1:50, 10, 5)); colnames(dat) <- LETTERS[1:5]; dat ## A B C D E ## 1 1 11 21 31 41 ## 2 2 12 22 32 42 ## 3 3 13 23 33 43 ## 4 4 14 24 34 44 ## 5 5 15 25 35 45 ## 6 6 16 26 36 46 ## 7 7 17 27 37 47 ## 8 8 18 28 38 48 ## 9 9 19 29 39 49 ## 10 10 20 30 40 50
比如，我们现在有个任务，把第2-4列的和大于75的行提取出来。首先，要对第2-4列逐行求总和，那么我们可以用apply()函数：
dat$jh234 <- apply(dat[ , 2:4], 1, sum) dat ## A B C D E jh234 ## 1 1 11 21 31 41 63 ## 2 2 12 22 32 42 66 ## 3 3 13 23 33 43 69 ## 4 4 14 24 34 44 72 ## 5 5 15 25 35 45 75 ## 6 6 16 26 36 46 78 ## 7 7 17 27 37 47 81 ## 8 8 18 28 38 48 84 ## 9 9 19 29 39 49 87 ## 10 10 20 30 40 50 90
可以看到这里增加了一列，这一列就是2-4列的和，apply()里第一个参数是要计算的列组成的子集，然后1表示逐行进行计算，sum表示求和。
然后我们在根据这一列来进行取子集。
Dat <- dat[dat$jh234 > 75, ] Dat ## A B C D E jh234 ## 6 6 16 26 36 46 78 ## 7 7 17 27 37 47 81 ## 8 8 18 28 38 48 84 ## 9 9 19 29 39 49 87 ## 10 10 20 30 40 50 90
同样的，我们可以求每一行的平均值(mean)、标准差(sd)等等的，比如这里对1-5列逐行求均值、标准差：
dat$平均 <- apply(dat[ , 1:5], 1, mean) dat$标准差 <- apply(dat[ , 1:5], 1, sd) dat ## A B C D E jh234 平均 标准差 ## 1 1 11 21 31 41 63 21 15.81139 ## 2 2 12 22 32 42 66 22 15.81139 ## 3 3 13 23 33 43 69 23 15.81139 ## 4 4 14 24 34 44 72 24 15.81139 ## 5 5 15 25 35 45 75 25 15.81139 ## 6 6 16 26 36 46 78 26 15.81139 ## 7 7 17 27 37 47 81 27 15.81139 ## 8 8 18 28 38 48 84 28 15.81139 ## 9 9 19 29 39 49 87 29 15.81139 ## 10 10 20 30 40 50 90 30 15.81139
另外，可以使用自定义的函数，先定义一个函数，然后用函数名，或者直接就写出函数，比如可以这样，自定义了一个求各个元素的平方的加和的函数：
pfh <- function(x) sum(x^2) dat$平方和 <- apply(dat[ , 1:5], 1, pfh) dat$平方和2 <- apply(dat[ , 1:5], 1, function(x) sum(x^2)) dat ## A B C D E jh234 平均 标准差 平方和 平方和2 ## 1 1 11 21 31 41 63 21 15.81139 3205 3205 ## 2 2 12 22 32 42 66 22 15.81139 3420 3420 ## 3 3 13 23 33 43 69 23 15.81139 3645 3645 ## 4 4 14 24 34 44 72 24 15.81139 3880 3880 ## 5 5 15 25 35 45 75 25 15.81139 4125 4125 ## 6 6 16 26 36 46 78 26 15.81139 4380 4380 ## 7 7 17 27 37 47 81 27 15.81139 4645 4645 ## 8 8 18 28 38 48 84 28 15.81139 4920 4920 ## 9 9 19 29 39 49 87 29 15.81139 5205 5205 ## 10 10 20 30 40 50 90 30 15.81139 5500 5500
使用apply可以逐列进行运算，这时候把第二个参数变成2即可。 比如这里要对1-5列逐列计算每列的平均值：
apply(dat[ , 1:5], 2, mean) ## A B C D E ## 5.5 15.5 25.5 35.5 45.5 apply(dat[ , 1:5], 2, function(x) sum(sqrt(x))) ## A B C D E ## 22.46828 39.19770 50.41687 59.53294 67.42001
lapply序列
lapply序列的函数包括lapply、sapply、vapply、replicate、simplify2array。
先构建一个列表用于本节的示例：
L <- list(A=1:3, B=4:6, C=1:6) L ## $A ## [1] 1 2 3 ## ## $B ## [1] 4 5 6 ## ## $C ## [1] 1 2 3 4 5 6
lapply的用法如下：
lapply(X, FUN, ...)
lapply对X的每个元素逐个执行函数FUN，返回一个同等长度的列表。X可以是向量或者列表。
lapply(L, sum) ## $A ## [1] 6 ## ## $B ## [1] 15 ## ## $C ## [1] 21
sapply
mapply
mapply是sapply的多变量版本，返回的值是向量, 注意这个mapply的使用方法，方程写在前边，这里的x，y，z可以用i,j,k之类的，可以有更多的参数。对于多个变量不是同一类型的时候，显然使用mapply就比apply更合适。
mapply(function(x,y,z) x + y^2 - z, x=1:10, y=11:20, z=21:30) ## [1] 101 124 149 176 205 236 269 304 341 380
tapply
tapply(X, INDEX, FUN = NULL, …, default = NA, simplify = TRUE)对对象X（通常是一个向量）根据INDEX分组，执行FUN指定的函数，…的位置提供FUN函数的更多的额参数。
x <- 1:10 set.seed(2023); (i <- sample(LETTERS[1:3], 10, replace=TRUE)) ## [1] "A" "C" "A" "C" "B" "B" "A" "A" "A" "A" tapply(X=x, INDEX=i, FUN=sum) ## A B C ## 38 11 6
by
by是tapply的一个面向对象的封装，用于数据框。用法为：by(data, INDICES, FUN, …, simplify = TRUE)。INDICES可以是向量，也可以是数据框。
CJ <- read.delim("ftp://ftp.biolab.wang/data/各科成绩.tsv", fileEncoding="GBK") by(data=CJ[ , c("科目", "成绩")], INDICES=CJ[ , "性别"], FUN=summary) ## CJ[, "性别"]: 男 ## 科目 成绩 ## Length:2200 Min. :31.00 ## Class :character 1st Qu.:61.00 ## Mode :character Median :67.00 ## Mean :67.47 ## 3rd Qu.:75.00 ## Max. :95.00 ## ------------------------------------------------------------------------------------------------------- ## CJ[, "性别"]: 女 ## 科目 成绩 ## Length:2648 Min. :36.00 ## Class :character 1st Qu.:63.00 ## Mode :character Median :70.00 ## Mean :69.79 ## 3rd Qu.:77.00 ## Max. :99.00 CJx <- CJ[CJ$科目 %in% c("细胞生物学", "遗传学"), ] by(data=CJx[ , c("成绩")], INDICES=CJx[ , c("科目", "性别")], FUN=summary) ## 目: 细胞生物学 ## 性别: 男 ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 40.00 61.00 66.00 67.15 73.00 87.00 ## ------------------------------------------------------------------------------------------------------- ## 科目: 遗传学 ## 性别: 男 ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 43.00 63.00 69.00 69.59 76.00 92.00 ## ------------------------------------------------------------------------------------------------------- ## 科目: 细胞生物学 ## 性别: 女 ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 44.00 62.50 69.00 69.37 76.00 89.00 ## ------------------------------------------------------------------------------------------------------- ## 科目: 遗传学 ## 性别: 女 ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 44.00 65.00 72.00 71.98 78.00 92.00
aggregate
aggregate(formula=y~x, FUN), 根据x将y分为多个组，每个组执行FUN。
aggregate(x, by, FUN, …, simplify = TRUE, drop = TRUE)将数据框x按照by分成不同的子集，然后执行FUN的函数。
x <- 1:10 set.seed(2023); (i <- sample(LETTERS[1:3], 10, replace=TRUE)) ## [1] "A" "C" "A" "C" "B" "B" "A" "A" "A" "A" aggregate(x~i, FUN=max) ## i x ## 1 A 10 ## 2 B 6 ## 3 C 4 aggregate(data.frame(x), by=data.frame(i), FUN=sum) ## i x ## 1 A 38 ## 2 B 11 ## 3 C 6
表格形状转换
我们先创建这个示例的数据框，宽型的：
dat_wide <- data.frame( 人员 = c("张三", "李四", "王五"), 年龄 = c(50, 60, 70), 舒张压 = c(80, 85, 90), 收缩压 = c(120, 135, 150), 心率 = c(72, 75, 68) ) dat_wide ## 人员 年龄 舒张压 收缩压 心率 ## 1 张三 50 80 120 72 ## 2 李四 60 85 135 75 ## 3 王五 70 90 150 68
使用reshape函数进行转换，有几个参数这里简要讲解，还有更复杂的用法用到了在去读手册。这里简单说能实现我们这种目的的用法：
1. data 参数是我们要转换的数据框，直接放在第一个就行，不用写参数名。后边为了方便我们更改了参数的顺序，就写出参数名了。
2. direction参数设置方向，宽型转长型设置direction = “long”，我比较喜欢把它放到第二个来，方便看出目的。
3. varing 指定哪些变量列是需要合并的，可以是列名，也可以是列的序号。比如这里可以写成varying = 3:5。
4. times 指定这些变量名称放到一列之后对应的值，比如我这里直接就用列名了，如有必要可以改，比较灵活。
5. timevar 长型数据框中，存储变量名称一列的列名，这一列用来区分某一行对应原始的哪一列。
6. idvar 指定标识来自原先同一行（观测对象）的列的名称，要求能够重复值，如果一列不够唯一就用更多的多列，不然会报错。如果并不设置这一参数，最后会多出一列叫id的整数型的列，用来标识来自同一行的新行。
dat_long <- reshape(dat_wide, direction = "long", varying = 3:5, times = names(dat_wide)[3:5], v.names = "数值", timevar = "指标", idvar = "人员") dat_long ## 人员 年龄 指标 数值 ## 张三.舒张压 张三 50 舒张压 80 ## 李四.舒张压 李四 60 舒张压 85 ## 王五.舒张压 王五 70 舒张压 90 ## 张三.收缩压 张三 50 收缩压 120 ## 李四.收缩压 李四 60 收缩压 135 ## 王五.收缩压 王五 70 收缩压 150 ## 张三.心率 张三 50 心率 72 ## 李四.心率 李四 60 心率 75 ## 王五.心率 王五 70 心率 68 # 得到的数据框还带有"reshapeLong"的附加属性。可以如下命令查看 str(dat_long) attributes(dat_long)
接下来我们在看看把这个数据框转回去成为款型的。这回没有新的参数需要介绍，只要改变方向，设置三个参数：
dat_wide_again <- reshape(dat_long, direction = "wide", timevar = "指标", v.names = "数值", idvar = "人员") dat_wide_again ## 人员 年龄 数值.舒张压 数值.收缩压 数值.心率 ## 张三.舒张压 张三 50 80 120 72 ## 李四.舒张压 李四 60 85 135 75 ## 王五.舒张压 王五 70 90 150 68 str(dat_wide_again) # 新的列名里边含有“数值.”这个前缀，如果想要去掉，可以这样 names(dat_wide_again) <- gsub("^数值[.]", "", names(dat_wide_again))
很多年前，reshape::melt()和reshape::cast()是比较先进的函数，我也比较喜欢使用，用起来也方便。现在我更倾向于教大家直接用这个stats::reshape()函数了，因为这个reshape()函数更底层一些，用起来更灵活。 如果反复要用，喜欢简单点，可以自己封装一个函数。比如这样子：
# 我们来封装一个函数，只要指定用作id的列，其余的列当做要合并的变量，这个函数相当于reshape::melt() W2L <- function(x, id){ X <- reshape(x, direction = "long", varying = names(x)[-id], times = names(x)[-id], v.names = "value", timevar = "variable") X$id <- NULL row.names(X) <- NULL data.frame(X) } W2L(dat_wide, id=1:2) ## 人员 年龄 variable value ## 1 张三 50 舒张压 80 ## 2 李四 60 舒张压 85 ## 3 王五 70 舒张压 90 ## 4 张三 50 收缩压 120 ## 5 李四 60 收缩压 135 ## 6 王五 70 收缩压 150 ## 7 张三 50 心率 72 ## 8 李四 60 心率 75 ## 9 王五 70 心率 68
# 我们来封装一个函数，只要指定用作id、variable、value列的名称，就能从长型的转为宽型的，这个函数相当于reshape::cast() L2W <- function(x, id, variable, value){ X <- reshape(x[ , c(id, variable, value)], direction = "wide", timevar = variable, v.names = value, idvar = id) row.names(X) <- NULL data.frame(X) } L2W(dat_long, id="人员", variable="指标", value="数值") ## 人员 数值.舒张压 数值.收缩压 数值.心率 ## 1 张三 80 120 72 ## 2 李四 85 135 75 ## 3 王五 90 150 68 L2W(dat_long, id=c("人员", "年龄"), variable="指标", value="数值") ## 人员 年龄 数值.舒张压 数值.收缩压 数值.心率 ## 1 张三 50 80 120 72 ## 2 李四 60 85 135 75 ## 3 王五 70 90 150 68
表格合并
先构建两个示例的数据框
A <- data.frame(X=LETTERS[1:5], Y=letters[1:5], Z=1:5) A ## X Y Z ## 1 A a 1 ## 2 B b 2 ## 3 C c 3 ## 4 D d 4 ## 5 E e 5 B <- data.frame(X=LETTERS[c(1, 3, 5, 6, 8)], Y=letters[6:10], Z=6:10) B ## X Y Z ## 1 A f 6 ## 2 C g 7 ## 3 E h 8 ## 4 F i 9 ## 5 H j 10
rbind
表格合并有多种，一种是A和B两个表头相同，将B增补到A的后边，这时候使用rbind函数。 rbind的时候，需要两个数据框的表头(names)一致。
C <- rbind(A, B) C ## X Y Z ## 1 A a 1 ## 2 B b 2 ## 3 C c 3 ## 4 D d 4 ## 5 E e 5 ## 6 A f 6 ## 7 C g 7 ## 8 E h 8 ## 9 F i 9 ## 10 H j 10
将两个表头相同的数据框合并，B增补到A的后边
D <- rbind(A, B, C) D ## X Y Z ## 1 A a 1 ## 2 B b 2 ## 3 C c 3 ## 4 D d 4 ## 5 E e 5 ## 6 A f 6 ## 7 C g 7 ## 8 E h 8 ## 9 F i 9 ## 10 H j 10 ## 11 A a 1 ## 12 B b 2 ## 13 C c 3 ## 14 D d 4 ## 15 E e 5 ## 16 A f 6 ## 17 C g 7 ## 18 E h 8 ## 19 F i 9 ## 20 H j 10
rbind可以同时将多个数据框rbind在一起，按顺序增补。
rbind可以用于任意数据框和0行0列的数据框，这个特性就很好，可以用于数据框的初始化，比如用于for循环中每次迭代rbind一个/些表格。
E <- data.frame() E ## data frame with 0 columns and 0 rows for(i in 1:3){ E <- rbind(E, A) } E ## X Y Z ## 1 A a 1 ## 2 B b 2 ## 3 C c 3 ## 4 D d 4 ## 5 E e 5 ## 6 A a 1 ## 7 B b 2 ## 8 C c 3 ## 9 D d 4 ## 10 E e 5 ## 11 A a 1 ## 12 B b 2 ## 13 C c 3 ## 14 D d 4 ## 15 E e 5
cbind
另外是将两个数据框并排放放置，这时候要求两个数据框的行数是一致的。cbind的时候，两个数据框的行名可以不一致，最后保留的是前边一个的行名。使用cbind的时候需要小心，要知道两个数据框行的排列顺序是否需要注意。
先将A和B的表头改一下，以便后边观察：
E <- A; names(E) <- c("X1", "Y1", "Z1") F <- A; names(F) <- c("X2", "Y2", "Z2")
G <- cbind(E, F) G ## X1 Y1 Z1 X2 Y2 Z2 ## 1 A a 1 A a 1 ## 2 B b 2 B b 2 ## 3 C c 3 C c 3 ## 4 D d 4 D d 4 ## 5 E e 5 E e 5 H <- cbind(E, F, G) H ## X1 Y1 Z1 X2 Y2 Z2 X1 Y1 Z1 X2 Y2 Z2 ## 1 A a 1 A a 1 A a 1 A a 1 ## 2 B b 2 B b 2 B b 2 B b 2 ## 3 C c 3 C c 3 C c 3 C c 3 ## 4 D d 4 D d 4 D d 4 D d 4 ## 5 E e 5 E e 5 E e 5 E e 5
merge
很多时候，我们需要将两个数据框按照特定的列来对应着合并，这时候就需要用到merge函数。比如我们有两个表格，一列表格记录学生的期考成绩，一列表格记录学生的平时成绩，那要把这两个成绩放在一个表格中，那么我们需要按照学号和姓名对应起来，不能简单的用cbind，因为两个表格里边学号的排列顺序可能不一致，有的学生可能在其中一个表格里还没有。
AB <- merge(A, B, by="X") AB ## X Y.x Z.x Y.y Z.y ## 1 A a 1 f 6 ## 2 C c 3 g 7 ## 3 E e 5 h 8
根据X列来merge数据框A和B，可以看到，只有X列两个表格都有的元素被保留了。同时也看到了，两个表格相同表头的列的名称被加上了.x和.y以示区别。
可以设置，这个后缀:
AB <- merge(A, B, by="X", suffixes = c(".A",".B")) AB ## X Y.A Z.A Y.B Z.B ## 1 A a 1 f 6 ## 2 C c 3 g 7 ## 3 E e 5 h 8
可以控制两个表格的保留情况：
AB <- merge(A, B, by="X", all.x=TRUE) AB ## X Y.x Z.x Y.y Z.y ## 1 A a 1 f 6 ## 2 B b 2 <NA> NA ## 3 C c 3 g 7 ## 4 D d 4 <NA> NA ## 5 E e 5 h 8 AB <- merge(A, B, by="X", all.y=TRUE) AB ## X Y.x Z.x Y.y Z.y ## 1 A a 1 f 6 ## 2 C c 3 g 7 ## 3 E e 5 h 8 ## 4 F <NA> NA i 9 ## 5 H <NA> NA j 10 AB <- merge(A, B, by="X", all=TRUE) AB ## X Y.x Z.x Y.y Z.y ## 1 A a 1 f 6 ## 2 B b 2 <NA> NA ## 3 C c 3 g 7 ## 4 D d 4 <NA> NA ## 5 E e 5 h 8 ## 6 F <NA> NA i 9 ## 7 H <NA> NA j 10
AE <- merge(A, E, by.x = c("X"), by.y = c("X1")) AE ## X Y Z Y1 Z1 ## 1 A a 1 a 1 ## 2 B b 2 b 2 ## 3 C c 3 c 3 ## 4 D d 4 d 4 ## 5 E e 5 e 5
当两个表格用于合并的一列的表头不一样的时候，可以用by.x和by.y分别设置。
表格分组统计
比如现在我们有一个如下一个任务，希望对几个指标分类别统计，如何来实现？
一个非常直观的做法就是把表格按照分裂拆分成多个子表格，然后每个表格进行计算，算好之后再把结果合并。
接下来我们用R来一步一步实现。首先我们先做一个示例数据框，并且给它转成长型的。
dat_wide <- data.frame( 姓名 = c("张三", "李四", "王五", "赵六", "钱七"), 性别 = factor(c("男", "女", "男", "女", "男")), 舒张压 = c(78, 72, 82, 75, 80), 收缩压 = c(120, 115, 128, 118, 125), 心率 = c(72, 68, 75, 70, 78) ) dat_wide ## 姓名 性别 舒张压 收缩压 心率 ## 1 张三 男 78 120 72 ## 2 李四 女 72 115 68 ## 3 王五 男 82 128 75 ## 4 赵六 女 75 118 70 ## 5 钱七 男 80 125 78 dat_long <- reshape(dat_wide, direction = "long", varying = c("舒张压", "收缩压", "心率"), v.names = "数值", timevar = "指标", times = c("舒张压", "收缩压", "心率"), idvar = "姓名" ) dat_long ## 姓名 性别 指标 数值 ## 张三.舒张压 张三 男 舒张压 78 ## 李四.舒张压 李四 女 舒张压 72 ## 王五.舒张压 王五 男 舒张压 82 ## 赵六.舒张压 赵六 女 舒张压 75 ## 钱七.舒张压 钱七 男 舒张压 80 ## 张三.收缩压 张三 男 收缩压 120 ## 李四.收缩压 李四 女 收缩压 115 ## 王五.收缩压 王五 男 收缩压 128 ## 赵六.收缩压 赵六 女 收缩压 118 ## 钱七.收缩压 钱七 男 收缩压 125 ## 张三.心率 张三 男 心率 72 ## 李四.心率 李四 女 心率 68 ## 王五.心率 王五 男 心率 75 ## 赵六.心率 赵六 女 心率 70 ## 钱七.心率 钱七 男 心率 78
接下来用split()函数将这个表格按照第2-3列拆分，这两列内容相同的行被拆到同一个子表格：
dat_list <- split(dat_long, dat_long[ , 2:3]) dat_list ## $男.收缩压 ## 姓名 性别 指标 数值 ## 张三.收缩压 张三 男 收缩压 120 ## 王五.收缩压 王五 男 收缩压 128 ## 钱七.收缩压 钱七 男 收缩压 125 ## ## $女.收缩压 ## 姓名 性别 指标 数值 ## 李四.收缩压 李四 女 收缩压 115 ## 赵六.收缩压 赵六 女 收缩压 118 ## ## $男.舒张压 ## 姓名 性别 指标 数值 ## 张三.舒张压 张三 男 舒张压 78 ## 王五.舒张压 王五 男 舒张压 82 ## 钱七.舒张压 钱七 男 舒张压 80 ## ## $女.舒张压 ## 姓名 性别 指标 数值 ## 李四.舒张压 李四 女 舒张压 72 ## 赵六.舒张压 赵六 女 舒张压 75 ## ## $男.心率 ## 姓名 性别 指标 数值 ## 张三.心率 张三 男 心率 72 ## 王五.心率 王五 男 心率 75 ## 钱七.心率 钱七 男 心率 78 ## ## $女.心率 ## 姓名 性别 指标 数值 ## 李四.心率 李四 女 心率 68 ## 赵六.心率 赵六 女 心率 70
接下来用lapply()函数逐个对这些表格进行处理：
res_list <- lapply(dat_list, function(x){ res_x <- data.frame( 平均值 = sum(x$数值), 标准差 = sd(x$数值), 个体数 = length(x$数值) ) res_x$均值标准误 = res_x$标准差/sqrt(res_x$个体数) dat_x <- x[1 , 2:3] return(cbind(dat_x, res_x)) }) res_list ## $男.收缩压 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 张三.收缩压 男 收缩压 373 4.041452 3 2.333333 ## ## $女.收缩压 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 李四.收缩压 女 收缩压 233 2.12132 2 1.5 ## ## $男.舒张压 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 张三.舒张压 男 舒张压 240 2 3 1.154701 ## ## $女.舒张压 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 李四.舒张压 女 舒张压 147 2.12132 2 1.5 ## ## $男.心率 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 张三.心率 男 心率 225 3 3 1.732051 ## ## $女.心率 ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 李四.心率 女 心率 138 1.414214 2 1
接下来用do.call()函数将这些结果合并成一个表格：
res <- do.call("rbind", res_list) res ## 性别 指标 平均值 标准差 个体数 均值标准误 ## 男.收缩压 男 收缩压 373 4.041452 3 2.333333 ## 女.收缩压 女 收缩压 233 2.121320 2 1.500000 ## 男.舒张压 男 舒张压 240 2.000000 3 1.154701 ## 女.舒张压 女 舒张压 147 2.121320 2 1.500000 ## 男.心率 男 心率 225 3.000000 3 1.732051 ## 女.心率 女 心率 138 1.414214 2 1.000000
这样的分组统计我以往常用plyr::ddply()函数，但plyr包目前已经不再先进，其继任者dplyr属于tidyverse方言与基础R不同，这里不增加学习负担，用已学习的内容解决我们的问题即可。
如果数据量很大，比如分组很多，如果最后do.call(“rbind”)比较慢，可以使用data.table::rbindlist()函数。
屏幕输入
readline()从屏幕输入。
X <- readline(prompt="请输入一个字母：")
LiveDays <- function() { X <- readline(prompt="输入您的生日：") Y <- Sys.Date() - as.Date(X) print(sprintf("您已经出生%d天了。", Y)) } LiveDays()
其它常用函数
stack与split
stack将元素为向量的列表转为数据框
L <- list(A=1:3, B=5:6, C=10) D <- stack(L) D ## values ind ## 1 1 A ## 2 2 A ## 3 3 A ## 4 5 B ## 5 6 B ## 6 10 C
而split可以把上方的过程进行逆转，根据一列将另一列拆分到列表：
E <- split(D$values, D$ind) E ## $A ## [1] 1 2 3 ## ## $B ## [1] 5 6 ## ## $C ## [1] 10
seq
seq产生等差数列，从第一个参数开始，以第三个参数为差，不大于第二个数。使用length参数，产生指定长度的等差数列，差通过计算得到。
seq(2, 13, 2) ## [1] 2 4 6 8 10 12 seq(2, 13, length=6) ## [1] 2.0 4.2 6.4 8.6 10.8 13.0
rep
rep用对对向量进行重复指定的次数，用each对每个元素依次进行重复。
rep(c(1,3,4), 3) ## [1] 1 3 4 1 3 4 1 3 4 rep(c(1,3,4), each=3) ## [1] 1 1 1 3 3 3 4 4 4
unique
combn
combn用于产生组合：
combn(1:5, 2) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] ## [1,] 1 1 1 1 2 2 2 3 3 4 ## [2,] 2 3 4 5 3 4 5 4 5 5 combn(1:5, 3) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] ## [1,] 1 1 1 1 1 1 2 2 2 3 ## [2,] 2 2 2 3 3 4 3 3 4 4 ## [3,] 3 4 5 4 5 5 4 5 5 5
stopifnot
stopifnot函数用于判断一个逻辑向量是否都是TRUE
x <- 1:10 stopifnot(x <= 10) stopifnot(x <= 5) ## 错误: 不是所有的x <= 5都是TRUE
all
which
which返回逻辑向量中TRUE的值的index
x <- 1:10 which(x<=8 & x >=5) ## [1] 5 6 7 8 which(x %% 2 == 0) ## [1] 2 4 6 8 10
with
with用于指定表达式的数据环境，用法为with(data, expr)，其中data可以是数据框或者列表，expr中可以直接使用其元素的名称(数据框中为列的名称)来调用数据。
rm(list=ls()) X <- data.frame(A=sample(1:10), B=sample(3:12)) t.test(A, B) ## Error in t.test(A, B) : 找不到对象'A' with(X, t.test(A, B)) ## Welch Two Sample t-test ## ## data: A and B ## t = -1.4771, df = 18, p-value = 0.1569 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -4.8446619 0.8446619 ## sample estimates: ## mean of x mean of y ## 5.5 7.5
cut
cut(x, breaks)将一个数值型向量x按照breaks指定的界限cut成不同的范围段转的因子。
set.seed(2002); x <- sample(1:100, 10) cut(x, c(0, 10, 50, 101)) ## [1] (10,50] (50,101] (50,101] (50,101] (10,50] (50,101] (50,101] (50,101] (10,50] (10,50] ## Levels: (0,10] (10,50] (50,101]
scale
scale对一个矩阵各列进行标准化，
set.seed(2023); x <- sample(1:100) x <- matrix(x, nrow=5) x ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] ## [1,] 80 65 98 4 38 96 68 32 42 37 28 61 27 59 18 53 35 77 25 62 ## [2,] 47 29 63 91 41 48 34 46 75 60 97 31 16 100 10 23 55 2 13 43 ## [3,] 72 49 3 24 30 78 40 90 36 19 15 50 20 69 99 54 85 8 58 84 ## [4,] 26 81 79 9 39 33 52 86 17 56 21 82 67 92 66 6 73 7 51 1 ## [5,] 44 5 45 57 93 70 83 64 95 76 12 22 88 87 11 94 71 74 14 89 scale(x) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] ## [1,] 1.1918945 0.643152 1.1141701 -0.9016595 -0.4017180 1.2480130 0.6247123 -1.26218376 -0.3498681 -0.5719597 -0.1863559 0.49351222 -0.5176642 -1.3238456 -0.5709464 0.2076411 -1.4908078 1.1322296 -0.3421982 0.1735871 ## [2,] -0.3093467 -0.562758 0.1489237 1.4754429 -0.2835657 -0.6843942 -1.0610193 -0.70298842 0.6997362 0.4720937 1.7619106 -0.76117987 -0.8606948 1.0992647 -0.7712785 -0.6822493 -0.4555246 -0.8243883 -0.9125286 -0.3583733 ## [3,] 0.8279573 0.107192 -1.5057844 -0.3551992 -0.7167910 0.5233603 -0.7635373 1.05448264 -0.5407052 -1.3890449 -0.5534207 0.03345846 -0.7359564 -0.7328431 1.4574158 0.2373041 1.0974002 -0.6678589 1.2262103 0.7895413 ## [4,] -1.2646820 1.179112 0.5901792 -0.7650445 -0.3623339 -1.2882715 -0.1685732 0.89471254 -1.1450229 0.2905192 -0.3840062 1.37179668 0.7297195 0.6264627 0.6310460 -1.1865206 0.4762303 -0.6939472 0.8935176 -1.5342859 ## [5,] -0.4458231 -1.366698 -0.3474887 0.5464603 1.7644087 0.2012924 1.3684175 0.01597701 1.3358600 1.1983917 -0.6381279 -1.13758749 1.3845959 0.3309614 -0.7462370 1.4238247 0.3727020 1.0539648 -0.8650010 0.9295308 ## attr(,"scaled:center") ## [1] 53.8 45.8 57.6 37.0 48.2 65.0 55.4 63.6 53.0 49.6 34.6 49.2 43.6 81.4 40.8 46.0 63.8 33.6 32.2 55.8 ## attr(,"scaled:scale") ## [1] 21.98181 29.85297 36.26017 36.59918 25.39094 24.83948 20.16928 25.03597 31.44042 22.02953 35.41610 23.91025 32.06712 16.92040 39.93370 33.71202 19.31839 38.33145 21.04044 35.71694
scale(1:5) ## [,1] ## [1,] -1.2649111 ## [2,] -0.6324555 ## [3,] 0.0000000 ## [4,] 0.6324555 ## [5,] 1.2649111 ## attr(,"scaled:center") ## [1] 3 ## attr(,"scaled:scale") ## [1] 1.581139
sample (随机取样)
随机取样用sample函数。 这里从从1-100中随机选出10个数：
sample(1:100, 10) ## [1] 61 35 24 16 39 34 44 54 33 6 sample(1:100, 10) ## [1] 59 97 26 31 16 98 61 9 13 50 sample(1:100, 10) ## [1] 93 75 83 13 31 46 77 60 10 94
每次执行的结果都不一样， 以上三次输出再再次运行时可能都不再会碰到。要保证可重复，需要在sample之前set.seed,
set.seed(1) sample(1:100, 10) ## [1] 68 39 1 34 87 43 14 82 59 51 set.seed(1) sample(1:100, 10) ## [1] 68 39 1 34 87 43 14 82 59 51 set.seed(2022) sample(1:100, 10) ## [1] 100 51 78 55 75 68 6 63 98 14 set.seed(2022) sample(1:100, 10) ## [1] 100 51 78 55 75 68 6 63 98 14 sample(1:100, 10) ## [1] 97 85 21 54 74 7 73 79 99 37
set.seed里的参数是一个整数，设置同样的 seed的情况下，每次sample输出的随机数是一样的。有时候记录这个seed值以便实验能够被别人重复出来。不先设置seed，得到的随机抽样又不同了，最后这次得到的结果可能不在能够重复得到。
set.seed(2022); sample(c("a", "b", "c", "d", "e"), 3) ## [1] "d" "c" "b" set.seed(2022); sample(c("a", "b", "c", "d", "e"), 6) ## Error in sample.int(length(x), size, replace, prob) : ## 'replace = FALSE'，因此不能取比总体要大的样本 set.seed(2022); sample(c("a", "b", "c", "d", "e"), 6, replace=TRUE) ## [1] "d" "c" "c" "d" "d" "c" set.seed(2022); sample(c("a", "b", "c", "d", "e"), 6, replace=TRUE, prob = c(5, 2, 1, 1, 1)) ## [1] "e" "b" "a" "b" "a" "b"
默认是不放回抽样，当抽样个数多于总体个数时，就会报错。设置replace=TRUE就是放回抽样了，prob 参数设置各元素被抽到的概率权重，这里a被抽到的概率最高。
table(replicate(2000, sample(c("a", "b", "c", "d", "e"), 6, replace=TRUE, prob = c(5, 4, 1, 1, 1)))) ## a b c d e ## 5054 3967 1004 972 1003
统计重复2000次抽到各元素的总数，可以看到比例和我们设置的prob接近，重复次数越多，就会越接近这个比例。。
set.seed(2022); sample(1:100) ## [1] 100 51 78 55 75 68 6 63 98 14 7 97 86 69 1 99 32 65 29 12 ## [21] 79 60 44 70 93 22 87 11 20 62 3 95 82 66 34 39 33 31 19 56 ## [41] 88 71 53 50 41 61 77 52 84 10 64 35 26 72 48 16 23 76 5 45 ## [61] 2 24 83 17 81 15 18 42 13 91 94 80 85 92 27 8 25 49 74 36 ## [81] 43 73 9 54 37 4 89 30 21 40 57 96 38 67 46 59 90 47 58 28 set.seed(2022); sample(c("a", "b", "c", "d", "e")) ## [1] "d" "c" "b" "a" "e" set.seed(2022); sample(100) ## [1] 100 51 78 55 75 68 6 63 98 14 7 97 86 69 1 99 32 65 29 12 ## [21] 79 60 44 70 93 22 87 11 20 62 3 95 82 66 34 39 33 31 19 56 ## [41] 88 71 53 50 41 61 77 52 84 10 64 35 26 72 48 16 23 76 5 45 ## [61] 2 24 83 17 81 15 18 42 13 91 94 80 85 92 27 8 25 49 74 36 ## [81] 43 73 9 54 37 4 89 30 21 40 57 96 38 67 46 59 90 47 58 28
不设置抽取数目的时候默认抽取元素总数的数目, 对于抽取的列表只给一个整数的时候，默认是从1到此整数中来抽样。
? sample
了解更多的随机抽样的知识
正则表达式
正则表达式是一个用于表达一系列字符串的表达式(模式)，通过模式指定的模式匹配可以把包含这些字符串的字符串提取出来进行替换等等的。比如2020级的同学大家的学号都是11位，以2020开头，那么要找出一个文件中所有的学号就可以用这个信息去表示为2020[0-9]{7}[^0-9],这样匹配2020开头且有11位的数字，第12位就不能是0-9中的了。
用? regex查看更详细的信息。我这里按照我的理解梳理了一下，主要是便于用到的时候快速查找和回顾。若发现其中存在理解不到位或者翻译不准确的地方，欢迎告诉我，这里先行谢过。R使用两种正则表达式：
- 扩展的正则表达式(Extended Regular Expression, EREs)，默认情况下使用
- 类Perl正则表达式(Perl-like regular expression)，当使用grep等等函数时设置perl=TRUE时使用
扩展的正则表达式
R正则表达式中[ \ ^ $ . | ? * + ( ) {这12个字符是保留作特殊含义的，其本身不能被自己匹配，而其余的字符都是被能被自身所匹配的：
[]表示内部的任意字符将被匹配，连续的数字和字母可以用-连写，比如[0-9]等同于[0123456789], [a-d]等同于[abcd], [W-Z]等同于[WXYZ]；除了\和^以外，其它的10个保留的字符放到[]里边都是其原本的含义。^表示字符串起始，比如^A表示匹配模式处于字符串的开始，且以字母A开始。如果放到中括号中，表示意义相反，比如上边提到的[^0-9]就表示没有匹配到数字；$表示字符串结束，比如A$表示匹配模式处于字符串的末尾，且以字母A结尾。前文中的学号的例子，如果要匹配到字符串最末端的学号，那么就要把这种特殊的情况考虑进去，可以写成2020[0-9]{7}[^0-9]|2020[0-9]{7}$；|表示或者，比如上边的|两边的2020[0-9]{7}[^0-9]和2020[0-9]{7}$任意一个被匹配都算；.匹配除换行符以外的任意字符；*表示前边的字符或者字符组将被配0次或者多次(>=0次)；？表示前边的字符或者字符组将被配0次或者1次(0-1次)；+表示前边的字符或者字符组将被配1次或者多次(>=1次)；()表示内部的字符串作为一个整体被匹配；{}及其内部的数字表示前边的字符被匹配的次数，比如上边的[0-9]{7}就表示7位数字，每一位上都可以是0-9中任何一个。{n}表示被匹配n次，{n,}表示被匹配至少n次，{n,m}表示被匹配n至m次；\\表示转义。如果我们要匹配上边提到的这些字符本身，就需要使用转义，比如\\?匹配到?本身，而要匹配到\本身，则需要用[\\]匹配到\本身。因为\在字符串中就是具有特殊含义的，它本身要出现也是以\\的形式出现。
常规的字符前边加上\则依赖处理器而具有特定的含义，目前处理器如下这些：
\a表示BEL，也就是发出“咚”这样的一声，使用cat(“\a”)你就能够听到，当然要打开声音才可以\e表示ESC, 这个在我的电脑上似乎翻译不了\f表示FF, 换页\n表示LF, 换行\r表示CR, 回车\t表示TAB, 制表符
另外，还有一些预定义的模式，可以比较方便地使用，放到中括号[]里使用，比如[:alpha:]和[:digit:]。其中有这些：
[[:alnum:]]:[[:alpha:]]和[[:digit:]]，等同于[0-9A-Za-z][[:alpha:]]:[[:lower:]]和[[:upper:]]，等同于[A-Za-z][[:upper:]]: 大写字母，等同于[A-Z][[:lower:]]: 小写字母，等同于[a-z][[:digit:]]: 数字，即[0-9][[:xdigit:]]: 十六进制数字，即[0123456789ABCDEFabcdef]，等同于[0-9A-Fa-f][[:blank:]]: 空白字符，即：space(空格)、tab(制表符)、以及可能包括其它如非换行空格这样的locale依赖的字符[[:cntrl:]]：Control characters. In ASCII, these characters have octal codes 000 through 037, and 177 (DEL). In another character set, these are the equivalent characters, if any. 【没用到过，懒得翻译了】[[:graph:]]: Graphical characters:[[:alnum:]]and[[:punct:]].[[:print:]]: Printable characters:[[:alnum:]],[[:punct:]]and space.[[:punct:]]: 标点符号，即：! “ # $ % & ' ( ) * + , - . / : ; < = > ? @ [ \ ] ^ _ ` { | } ~[[:space:]]: 间隔字符，即制表符(tab,\t), 换行(newline,\n), 垂直制表符(vertical tab,\v), 换页符(form feed,\f), 回车(carriage return,\r), space and possibly other locale-dependent characters.
特殊含义字符的匹配
要匹配到方括号[, 可以把它放到方括号里，或者适用\\转义：
grep("[[]", "a[") ## [1] 1 grep("\\[", "a[") ## [1] 1
匹配\本身则需要放到方括号里并且适用双斜杠\\,字符串中的\也需要用双斜杠转义。
grep("[\\]", "a[\\") ## [1] 1 cat("A\\") ## A\
匹配^本身则需要用双斜杠转义：
grep("\\^", "a^") ## [1] 1
匹配$、.、|、?、*、+、(、)、{等特殊字符本身，以$为例，其余都一样：
grep("\\$", "a[\\^") ## integer(0) grep("[$]", "a[\\^") ## integer(0) grep("\\$", "a[\\^$") ## [1] 1 grep("[$]", "a[\\^$") ## [1] 1
预定义的模式的使用
比如，这里删掉(替换为无字符，也就是引号里什么也没有)所有的间隔字符或者所有的字母和数字：
gsub("[[:space:]]", "", "A B\tC\n1 2 3") ## [1] "ABC123" gsub("[[:alnum:]]", "", "A B\tC\n1 2 3") ## [1] " \t\n "
查看更多
? regex
R包安装
R之便利与强大，得益于大量的R包，通过安装新的R包获得新的功能。不同的R包可能来源不同，不同的来源，有不同的安装方法。这部分同学们不用都去尝试一遍，可以等自己需要安装包的时候再来参考。 至于哪个包来源于何处，一般上网搜索就能找到。有的包有很多来源，选择其一即可。一般来说，CRAN里边有的，就用CRAN的就可以了，安装也比较方便。
从存储库或本地文件安装(Install Packages from Repositories or Local Files)
这时候使用的是install.packages函数，使用? install.packages查看函数的用法、可以改变的参数以及更多的说明。其中常用的参数有pkgs、repos、lib、type，依次对应安装哪些包、从哪里获得安装文件、安装到哪里、安装文件是二进制包还是源代码。最简单的用法是指指定按住的包的名称，比如安装plyr包：
library(plyr) ## Error in library(plyr) : 不存在叫‘plyr’这个名字的程辑包 install.packages("plyr") ## --- 在此連線階段时请选用CRAN的鏡子 --- ## Warning: failed to download mirrors file (无法打开URL'https://cran.r-project.org/CRAN_mirrors.csv'); using local file 'D:/Program/R/R-4.3.1/doc/CRAN_mirrors.csv' ## 还安装相依关系‘Rcpp’ ## ## 试开URL’https://cloud.r-project.org/bin/windows/contrib/4.3/Rcpp_1.0.11.zip' ## Content type 'application/zip' length 2873801 bytes (2.7 MB) ## downloaded 2.7 MB ## ## 试开URL’https://cloud.r-project.org/bin/windows/contrib/4.3/plyr_1.8.8.zip' ## Content type 'application/zip' length 1161510 bytes (1.1 MB) ## downloaded 1.1 MB ## ## 程序包‘Rcpp’打开成功，MD5和检查也通过 ## 程序包‘plyr’打开成功，MD5和检查也通过 ## ## 下载的二进制程序包在 ## C:\Users\wanghong\AppData\Local\Temp\RtmpgnWGrl\downloaded_packages里 ## Warning message: ## In download.file(url, destfile = f, quiet = TRUE) : ## URL 'https://cran.r-project.org/CRAN_mirrors.csv': status was 'SSL connect error' library(plyr)
当一个包没有安装时候，如果使用library加载，就会报错，提示不存在这个包。安装成功之后，再library这个包就不会有报错信息了。如果正在使用的R比最新的R版本低的时候，新安装的R包可能是用更新版本的R来构建的，加载的时候会有提示，提示为警告(Warning)信息，一般不影响使用，一般可以忽略。
安装的时候，可以看到这些信息：
- 每次打开的R控制台第一次使用
install.packages安装时会弹出窗口让选择镜像，可以任选一个。如无法连接或太慢，可chooseCRANmirror()打开选择窗口重新选择镜像，可以选一个国内镜像，这样速度可能会快一些。另外，可以一次指定安装多个包，比如这里安装reshape、plyr、pheatmap：install.packages(c(“reshape”, “plyr”, “pheatmap”))。 - 有个警告
Warning: failed to download mirrors file …，说的是无法链接到r-project.org网页，这是网络问题，有时会遇到。安装R-base的时候已经附带了CRAN_mirrors.csv文件，这里就可以使用，所以不重新下载也能从这个文件中找到一个能够链接到的镜像。 - 接下来是限时安装依赖的包的名称，这些包也会被安装
- 接下来是下载这些包的提示，有时候比如因为网络不好之类的原因downloaded后边的大小不够，那么后边会安装不成功，再次运行安装命令即可。
- 接下会有
下载的二进制程序包在…里的提示，这里说的是安装包存放的位置，不是这个包安装的位置，这是临时文件，不用管它 - 接下来的警告和一开始的警告是一样的原因，并不是每次都会出现，这并不影响
为了避免选择镜像的麻烦，在安装的时候可以直接设置一个国内的镜像，可以用getCRANmirrors()函数获取镜像，选择.cn即可，比如这里从北大镜像下载：
install.packages("ggplot2", repos="https://mirrors.pku.edu.cn/CRAN/") ## 试开URL’https://mirrors.pku.edu.cn/CRAN/bin/windows/contrib/4.4/ggplot2_3.5.1.zip' ## Content type 'application/zip' length 5012242 bytes (4.8 MB) ## downloaded 4.8 MB ## ## 程序包‘ggplot2’打开成功，MD5和检查也通过 ## ## 下载的二进制程序包在 ## C:\Users\wanghong\AppData\Local\Temp\Rtmpqgr1Jl\downloaded_packages里
如上所示，不再提示让选择镜像了。
对于安装在哪里的问题，.Library是一个字符串，给定默认的程序库的路径。可以使用.libPaths函数查看当前的所有的包的路径：
.Library ## [1] "D:/Program/R/R-4.2.0/library" .libPaths() ## [1] "D:/Program/R/R-4.2.0/library"
默认情况下，以上两条命令应该得到同样的结果，此刻我的是D:/Program/R/R-4.2.0/library。如果要新增程序库的位置，可以使用.libPaths函数指定：
dir.create("D:/Program/R/R-4.2.0/library_test") .libPaths("D:/Program/R/R-4.2.0/library_test") .libPaths() ## [1] "D:/Program/R/R-4.2.0/library_test" "D:/Program/R/R-4.2.0/library"
这样，就能看到有两个路径了，而且新添加的是放在前边的。
如果不指定lib参数，那么在用install.packages安装R包的时候，默认就是安装到.libPaths()[1]对应的位置。有时候可能由于权限的问题，比如在集群上使用，自己没有权限在默认的库路径下安装包，这时候需要自己指定，就可以用这个.libPaths增加一个自己可读写的路径。这个路径可以直接用lib参数指定，但是在使用的时候依然需要用.libPaths来增加这个路劲，不然R是找不到这个包的。如下这两种安装路径的指定方式是一样的效果：
.libPaths("D:/Program/R/R-4.2.0/library_test") install.packages("ggplot2", repos="https://mirrors.pku.edu.cn/CRAN/")
.libPaths("") install.packages("ggplot2", lib=D:/Program/R/R-4.2.0/library_test, repos="https://mirrors.pku.edu.cn/CRAN/")
另外，就是type参数，有时候需要从源码安装，这时候就需要设定type=“source”，要能够从source安装，还需要安装Rtools。这个一般用不到，这里不展开了。需要用到的时候已经是有相当的水平了，自己也就能够弄明白了。
从Github安装
Github上的包的安装使用devtools包，用到install_github函数，比如这里安装github中的ggolot2。
if (!requireNamespace("devtools")) install.packages("devtools") devtools::install_github("tidyverse/ggplot2") #安装github中的开发中版本
一般在项目的页面上也会介绍如何安装的，根据介绍安装即可。
从Bioconductor安装
Bioconductor提供了许多高通量基因组数据的分析工具，对于来源于Bioconductor包的安装可以看具体页面的说明。比如用于基因差异分析的DESeq2网页上能够找到安装的方法：
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2")
如上，安装Bioconductor上的包用BiocManager包中的install函数。 如果没有安装过BiocManager包，则用前边提到的方法用install.packages函数把BiocManager包按上，详见官方安装指南。
如果网速比较慢，可以试试更改为国内的镜像，在安装之前先设置BioC_mirror选项，比如使用中国科技大学的镜像：
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") BiocManager::install("edgeR")
另外还能选择：
# 清华大学的镜像 https://mirrors.tuna.tsinghua.edu.cn/bioconductor/ # 南京大学的镜像 https://mirrors.nju.edu.cn/bioconductor/
有时候可能会遇到如下这样的警告信息，无法安装一些包，这时候可能是Bioconductor的版本过期了：
Warning message: packages ‘ComplexHeatmap’, ‘clusterProfiler’, ‘biomaRt’ are not available for Bioconductor version '3.17' Versions of these packages for your version of R might be available elsewhere, see the ideas at https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
可以这样更新到最新的版本：
BiocManager::install(version = "3.20")
从本地文件或者网络上的文件安装
通常用到的R包都会在以上的这些地方找到，有的可能不止一个地方存在，因此多种方式可以安装。极少数的情况，有的包并没有放到上边提到的这个软件库中，一般都会说明如何安装，比如把安装包下载下来按照。
有时候可能由于网络连接等问题，需要先把安装包或者源代码下载下来在安装，这时候也是使用install.packages函数。只是把pkgs参数设定成具体的安装包文件的路径、加上repos = NULL参数即可。安装包文件路径可以是网络链接。其余参数保持默认或者根据具体情况指定即可：
install.packages(“路径到文件或者网络连接到文件”, repos = NULL)
Github上的R包有时候或许是因为网络连接的问题，总是安装不成功，这时候可以下载到本地再安装。 当没有打包好的*.tar.gz这样的文件可供选择的时候，可以从Code下拉菜单Download ZIP，然后得到的XXX-master.zip这样的文件（XXX为包名称），这种情况下需要先把zip解压，解压之后得到XXX-master文件夹，文件夹重命名为XXX，把这个文件夹的路径作为参数，设置repos=NULL，type=“source”就可以了。
另外，可以参考这个InstaPrism包的安装示例。
从源码安装R程序包(需要安装Rtools)
从源码安装R包需要用到Rtools软件。Rtools不能通过常规install.packages()命令进行安装，需要到Rtools网页下载（根据自己R的版本选择合适的版本），像安装其他程序一样安装上。我是安装在D:\Program\R\rtools44文件夹的。安装好了之后，新打开一个R控制台，输入命令Sys.which(“make”)看看是否能够找到，如果没有则在R控制台需要执行如下代码，把安装路径添加到~/.Renviron文件中：
writeLines('PATH="D:\\Program\\R\\rtools44\\usr\\bin;${PATH}"', con = "~/.Renviron")
然后在R控制台执行如下命令试试看能否找到make的位置：
Sys.which("make") ## make ## "D:\\Program\\R\\rtools43\\usr\\bin\\make.exe"
这个命令应该能够看到make.exe文件的路径，如果看不到，那么说明没有配置好。可以再打开一个R Console或者重启一下R Console试试。
安装好了之后，就可以从源码安装R包了，也就是install.packages函数里再加上type=“source”选项。
R包的运用
R包的加载与使用
这里我们用ggplot2包为例演示R包的加载和探讨其中的一些事项。为了确保能够看到后边的演示效果，首先执行如下两行命令，确保ggplot2成功安装以及ggplot2没有预先加到搜索列表（search()命令查看）：
if(!requireNamespace("ggplot2")) install.packages("ggplot2", repos="https://mirrors.pku.edu.cn/CRAN/") if("package:ggplot2" %in% search()) detach(package:ggplot2, unload = TRUE)
要想使用R包里的函数，首先要把R包载入(load)到内存，要想能够直接通过函数的名称找到R包的函数，需要把R包附加（attach）到搜索列表。
library函数同时完成load和attach两个动作，加载R包的正确做法是使用library。 加载R包X使用library(X)命令，这样把R包X加载到内存，同时还把R包X的名称放到当前的工作空间的搜索列表（search()可查看）里，这样在调用X里边的函数的时候，就是从搜索列表里的R包里去找。后边加载的R包在搜索列表里边被放在前边，也就是如果有同名的函数，那么默认是用后加载的R包里的，因此需要注意。
search() ## [1] ".GlobalEnv" "package:stats" "package:graphics" "package:grDevices" "package:utils" ## [6] "package:datasets" "package:methods" "Autoloads" "package:base" # 没有加载ggplot2包，这时候在搜索列表里是没有没有package:ggplot2的 ggplot() + geom_point(aes(x=1, y=1)) ## Error in ggplot() : 没有"ggplot"这个函数 # 这时候找不到ggplot这个函数 library(gglot2) search() ## [1] ".GlobalEnv" "package:ggplot2" "package:stats" "package:graphics" "package:grDevices" ## [6] "package:utils" "package:datasets" "package:methods" "Autoloads" "package:base" # library载入R包之后，ggplot2就出现在了搜索列表中 ggplot() + geom_point(aes(x=1, y=1)) # 可以看到作图的窗口弹出来了，这时候就能直接用函数名调用ggplot2包里的函数了
require函数也能加载R包，但不要用它来加载R包。 require实际上就是在library函数外边加了个壳，使得在加载不成功的时候不像library那样报错，而是警告。加载成功则返回不可见的TRUE 同时返回FALSE（不可见，加载成功则）。因此，除非你确定需要用到它的返回值，要避免使用require来加载R包。
(library(ggplot2)) ## [1] "ggplot2" "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" (require(ggplot2)) ## [1] TRUE (library(ggplot3)) ## Error in library(ggplot3) : 不存在叫‘ggplot3’这个名字的程辑包 (require(ggplot3)) ## 载入需要的程辑包：ggplot3 ## [1] FALSE ## Warning message: ## In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : ## 不存在叫‘ggplot3’这个名字的程辑包
loadNamespace函数将R包载入内存而不附加到搜索列表，load而不attach。仅仅为了检查是否已经成功安装了某个R包，则可以使用requireNamespace函数，成功则返回TRUE，而不成功则是FALSE。这正如上边我们在安装之前有个if(!requireNamespace(“X”))这样的判断。
如果不把R包attach到搜索列表，调用R包里的函数就需要用ggplot2::ggplot()这样命令，在在函数名前边加上包的名称，中间用两个半角的冒号::连接，这样就能找到ggplot这个函数，而且明确使用的就是ggplot2里的ggplot这个函数，而不是别的包里的同名函数。
if("package:ggplot2" %in% search()) detach(package:ggplot2, unload = TRUE) # 为了保证效果，如果已经attach了则detach一下 ggplot() + geom_point(aes(x=1, y=1)) ## Error in ggplot() : 没有"ggplot"这个函数 ggplot2::ggplot() + geom_point(aes(x=1, y=1)) ## Error in geom_point(aes(x = 1, y = 1)) : 没有"geom_point"这个函数 ggplot2::ggplot() + ggplot2::geom_point(aes(x=1, y=1)) ## Error in aes(x = 1, y = 1) : 没有"aes"这个函数 # 没有把ggplot2加到搜索列表，就找不到ggplot2包里的ggplot、geom_point、aes等函数，加命名空间ggplot2之后就能找到了 ggplot2::ggplot()+ggplot2::geom_point(ggplot2::aes(x=1, y=1)) # 这时候就看到图出来了
如果一开始R包没有载入内存，那么使用ggplot2::ggplot()这种方式使R包被加载，执行一次之后在sessionInfo()的结果中能够看到这个包，列在loaded via a namespace (and not attached):栏目下，load了，但没attach，下次调用的时候，依然要这样在前边加上包的名称。在library(ggplot2)之后，我们会看到sessionInfo()结果里边多了许多列在loaded via a namespace (and not attached):栏目下的包，也就是依赖的包也会被载入内存，但是并没有附加到搜索列表，使用search()是看不到那些包的。
最新library的R包会被加到搜索列表的前边，也就是在搜索函数的时候优先从最新加载的R包里边找。因此，当不同的包里边含有相同名称的函数时，后library的R包里的函数就会掩盖（mask）先前载入的R包里的同名函数。比如reshape2包里边的melt等函数和reshape里的同名，library时就会有提示：
if(!requireNamespace("reshape")) install.packages("reshape", repos="https://mirrors.pku.edu.cn/CRAN/") if(!requireNamespace("reshape2")) install.packages("reshape2", repos="https://mirrors.pku.edu.cn/CRAN/") if("package:reshape" %in% search()) detach(package:reshape, unload = TRUE) if("package:reshape2" %in% search()) detach(package:reshape2, unload = TRUE) library(reshape) library(reshape2) ## ## 载入程辑包：‘reshape2’ ## ## The following objects are masked from ‘package:reshape’: ## ## colsplit, melt, recast library(reshape) search() ## [1] ".GlobalEnv" "package:reshape2" "package:reshape" "package:stats" "package:graphics" "package:grDevices" ## [7] "package:utils" "package:datasets" "package:methods" "Autoloads" "package:base"
很显然，先library了reshape，后library了reshape2，reshape2就放在了前边，即便再次执行library(reshape)并不会改变这个顺序了。
这时候就要注意，要知道你用的是哪个包的函数，不然就会弄错。这时候可以reshape::melt()、reshape2::melt()这样加上命名空间准确调用。为了避免对基础函数的掩盖，在使用与基础包里的函数名有冲突的包的时候，可以使用loadNamespace载入内存，然后加上包名称引用，以避免造成混淆。比如像被有的人吹捧的dplyr包就有些函数和基础包里的函数同名，这种冲突显然可能是作者欠考虑了。
可使用package ? X这样的命令查看R包X的使用指南，打开其说明文档，点页面下边的index可以看到目录，包括所有函数的页面链接：
package ? ggplot2 help(package="ggplot") # 这样效果是一样的
要了解一个包更多的信息，可以使用browseVignettes：
browseVignettes("ggplot2") # 查看手册
查看包是否已经加载，则用search函数：
search() # 查看已经加载的R包以及对象
要将一个R包移出搜索列表则，则使用detach函数：
detach("package:pheatmap", unload=TRUE) # 卸载已加载的R包，在需要换一个版本的包的时候可以用到
接下来会对一些包当中的常用函数做介绍，这些包的功能比这里介绍的要强大地多，根据自己的需要可以做出更多的学习。
R包运用示例之Reshape
这里原来的内容被移到了reshape，我现在推荐使用R核心包里stats的reshape函数，详见教程表格形状转换部分。stats::reshape函数更灵活，也足够方便，而且不用担心什么时候这个包就不再继续更新不兼容了。
R包运用示例之plyr
R包运用示例之 stringr
使用stringr包来比较方便地进行文本处理，安装加载后即可使用。
install.packages('stringr') library(stringr)
R包运用示例之 pheatmap
pheatmap是我比较常用得一个绘画热图得软件，比如拿到学生一个学期得课程成绩，那么我想看看这些课程成绩之间得相关性如何，我就可以使用热图展示出来。
library(pheatmap) library(reshape) options(scipen=20) dat <- read.delim("ftp://ftp.biolab.wang/data/成绩.tsv", header=TRUE, sep='\t', encoding="GBK") str(dat) ## 'data.frame': 830 obs. of 7 variables: ## $ 学号 : num 71024000001 71024000001 71024000001 71024000001 71024000002 ... ## $ 姓名 : chr "梁全" "梁全" "梁全" "梁全" ... ## $ 班级 : chr "3班" "3班" "3班" "3班" ... ## $ 课程编号: int 1006001 1006017 1006036 1006002 1006001 1006017 1006036 1006002 1006001 1006017 ... ## $ 课程名称: chr "体育" "解剖学" "物理学" "英语" ... ## $ 学分 : num 2 4.5 4 4 2 4.5 4 4 2 4.5 ... ## $ 成绩 : int 76 81 61 61 77 83 75 67 80 72 ... Dat <- data.frame(cast(dat, 学号+姓名~课程名称, value="成绩", fun.aggregate=max, fill=0)) str(Dat) ## 'data.frame': 161 obs. of 8 variables: ## $ 学号 : num 71024000001 71024000002 71024000003 71024000004 71024000005 ... ## $ 姓名 : chr "梁全" "李海丽" "李昕婧" "廖娇" ... ## $ 基础化学 : num 0 0 0 0 0 0 0 0 0 0 ... ## $ 解剖学 : num 81 83 72 80 77 81 68 81 84 82 ... ## $ 体育 : num 76 77 80 82 79 86 84 82 84 87 ... ## $ 物理学 : num 61 75 68 78 73 81 68 79 79 76 ... ## $ 细胞生物学: num 0 0 0 0 0 0 0 0 0 0 ... ## $ 英语 : num 61 67 81 83 76 70 81 70 83 83 ... COR <- cor(Dat[ , -c(1,2)], method= "spearman") COR ## 基础化学 解剖学 体育 物理学 细胞生物学 英语 ## 基础化学 1.0000000 0.13560152 0.13514623 -0.17063920 0.87764141 0.4830785 ## 解剖学 0.1356015 1.00000000 0.17643968 0.11719268 0.07082217 0.3650303 ## 体育 0.1351462 0.17643968 1.00000000 -0.08994174 0.15617101 0.1747764 ## 物理学 -0.1706392 0.11719268 -0.08994174 1.00000000 -0.21263873 -0.4220190 ## 细胞生物学 0.8776414 0.07082217 0.15617101 -0.21263873 1.00000000 0.4556115 ## 英语 0.4830785 0.36503030 0.17477638 -0.42201896 0.45561152 1.0000000 class(COR)
这时候就可以用heatmap 来看各个课程之间学生得分的相关性了
pheatmap(COR, main="Spearman Correlation")
很容易就看出物理学和英语竟然是负相关的， 而基础化学和细胞生物学的成绩相关性很高
COR <- cor(Dat[ , -c(1,2)], method= "pearson") pheatmap(COR, main="Pearson Correlation")
COR <- cor(Dat[ , -c(1,2)], method= "kendall") pheatmap(COR, main="Kendall Correlation")
pheatmap还有更为强大及丰富的功能，可查看使用手册
? pheatmap
R包运用示例之 ggplot2
ggplot2是很强大的作图R包，加上其强大的拓展, 学了ggplot2之后似乎都可以把基本的画图函数忘掉了。
以下是一些例子，可以试试结合各种元素画出不同风格和外观的图。
library(ggplot2) library(reshape) library(plyr) library(RColorBrewer) # 这里直接读入网络上的文件 my_expr_matrix <- read.delim("ftp://ftp.biolab.wang/data/rpkm.tsv", header=TRUE, sep='\t') # 这里的rpkm.tsv为归一化后的基因表达数据。 dat <- melt(my_expr_matrix, id=1) #用 ? melt 查看melt的用法 names(dat)[2:3] <- c("sample", "rpkm") dat$group <- gsub("[0-9]$", "", dat$sample) #样本去掉后边地数字是分组名（示例用，随便选的） dat$group <- factor(dat$group, levels=c("C","Tx","Ty")) #改变group的顺序 g <- #将后边的ggplot返回的对象“命名”为g ggplot() + geom_density(data=dat, aes(x=log2(rpkm+0.01), colour=sample)) g #使用g来喊出图来 #可以继续对g进行修饰、增加图层等等的 g + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) + theme_classic() # 从density图上可以看出log2(rpkm+0.01)=-5差不多可以作为基因是否表达的一个阈值。 # 可以直接画图，不存 ggplot() + geom_boxplot(data=dat, aes(x=sample, y=log2(rpkm+0.01))) # 如果把不表达的去掉，比如，至少需要有一个组中的至少两个样本表达超过这个阈值 Ns <- ddply(dat, .(geneID, group), summarize, Ns = sum(rpkm > 2^(-5) - 0.01)) # 求表达值超过阈值的样本数 Ng <- ddply(Ns, .(geneID), summarize, Ng = sum(Ns >= 2)) # 求至少两个样本超阈值的组数 Dat <- dat[dat$geneID %in% Ng[Ng$Ng > 0, ]$geneID, ] #仅保留有表达的基因 ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log2(rpkm+0.01))) #看看有什么不一样了 ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log2(rpkm+0.01), colour=sample)) # 增加colour ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log2(rpkm+0.01), colour=sample, fill=group)) #增加fill ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log2(rpkm+0.01), colour=group)) + scale_color_manual(values=c("red","green", "blue")) # 手动设置颜色 # 通过factor来改变样本展示的顺序 Dat$sample <- factor(Dat$sample, levels=c("Tx1", "Ty3", "Tx3", "C3", "C1", "Tx2", "C2", "Ty2", "Ty1")) ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log2(rpkm+0.01), colour=group)) + scale_color_manual(values=c("red","green", "blue")) ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log(rpkm+0.01), colour=group)) + scale_fill_manual(values=c("red","green", "blue")) + scale_x_discrete(labels=c(1:9)) #手动更改设置横坐标标记 g <- ggplot() + geom_boxplot(data=Dat, aes(x=sample, y=log(rpkm+0.0001), colour=group))+ xlab(NULL)+ ylab("log2(RPKM+0.0001)") + ggtitle("Boxplot of gene expression") + scale_y_continuous(limits=c(-10,10), breaks=seq(-10,10,by=1), expand=c(0,0))+ scale_colour_manual(values=c("red","green", "blue")) + theme_bw() + theme(legend.position="none", # ? theme查看更多的theme参数 plot.title=element_text(hjust=0.5, size=25), axis.text.x=element_text(size=15) ) ggsave("boxplot.pdf", plot=g) g <- ggplot() + geom_histogram(data=Dat, aes(x=log(rpkm+0.0001), fill=sample), position='dodge',binwidth=1) + #控制bin的宽度 facet_grid(.~group) + xlab("log2(RPKM+0.0001)")+ ggtitle("Histogram of gene expression") + scale_fill_manual(values=brewer.pal(9, 'Set1'), name='sample') + #用 ? RColorBrewer 查看这个颜色的信息 theme_bw() + theme(legend.position="none", # ? theme查看更多的theme参数 plot.title=element_text(hjust=0.5, size=25), axis.text.x=element_text(size=15) ) ggsave("histogram.pdf", plot=g) # 看一下两个组的表达均值的一致性 Data <- ddply(Dat, .(geneID, group), summarize, mean_rpkm=mean(rpkm)) Data <- cast(Data, geneID~group, value="mean_rpkm") pdf("散点图.pdf", family="GB1") #保持图片的时候中文可能会无法显示，使用family="GB1"即可解决这个问题 ggplot() + geom_point(data=Data, aes(x=log2(C+0.0001), y=log2(Tx+0.0001))) + xlab("对照")+ ylab("处理 X") + theme_bw() dev.off() # 散点图中包括圆圈(shape=1,10,13,16,19,20等)时，在一些可能阅读器中可能会显示为q，这时设置, 在ggsave中：useDingbats=FALSE \\
更详细的信息参见作图
在R里运行其它编程语言代码
Rcpp
使用Rcpp可能很容易的在R里边使用C++的代码，对一些需要优化执行效率的，可以使用C来写，然后借助Rcpp来调用，提高效率。需要安装Rcpp和Rtools。Rtools的安装见上文，而Rcpp的安装则和其它的大多数包一样：
install.packages("Rcpp", repos=getCRANmirrors()$URL[1])
同样是上文中的水仙花数的计算，我们可以用C来写，借助Rccp和Rtools进行编译之后在R里边运行。如下：
library(Rcpp) cppFunction('int SXHS(int n){ int N=1; for(int i=1; i<n; i++){ N *= 10; } for(int j=N; j<N*10; j++){ int S=0; int I=j; int s; for(int k=1; k<=n; k++){ s = I%10; int ss=1; for(int l=1; l<=n; l++){ ss *= s; } S += ss; I /= 10; } if(S==j){ Rprintf("%d \\n", j); } } return 0; }') SXHS(3) system.time(SXHS(8))
用C来算，像算N次方这样的都需要使用循环，而对于R则已经有现成的运算符，即便是这样，C依然是要快许多。这里使用R语言，算8位用了356s，而同样的电脑上，用C来写，算9位也只用了31s。
然而在大多数场合，我们不需要那么快，R的速度足够满足要求，而写代码的时间却大大节约了。写代码的时间才是我们的生命，而电脑运行的时间我们可以干别的。 而且，这里大家也看到了，通过Rcpp，能够很方便地把一段C代码整合到R里边。 对于一些关键的影响速度的环节，可以找个C高手改写一下代码。
JuliaConnectoR
JuliaConnectoR提供R和Julia之间的界面，从而可以从R运行Julia的脚本并且在两个语言间调用对象。Julia是近些年发展起来的一个面向科学计算的高性能动态语言，像R语言一样易学，而性能可与C这样的编译语言媲美。安装JuliaConnectoR之前先要安装Julia,比如我现在安装的是1.5.3版，安装在D:/Program/Julia/Julia-1.5.3。
install.packages("JuliaConnectoR", repos=getCRANmirrors()$URL[1])
library(JuliaConnectoR) juliaSetupOk() # 检查一下Julia设置是否OK，如果这里结果是TRUE，那么就可以用了。 juliaCall("sin.", c(1,2,3)) # 用Julia来计算1，2，3这几个数的正弦值
比如上文提到的水仙花数的寻找，我们可以使用Julia语言来写，然后从R里边调用：
Sxhs <- juliaEval(' function Sxhs!(n::Integer) for j::Integer in (10^(n-1)):(10^n - 1) S::Integer=0 I::Integer=j for k in 1:n S=S + (I%10)^n I=I ÷ 10 end if(S==j) println(S) end end end ') Sxhs(8)
这里发现从R里边运行的速度不如Julia里边运行那么快。我现在对Julia还不是很熟悉，所以无法给出解释。但是很显然，Julia会比R快一些。如果是在Julia的控制台运行，那速度会快很多。如下：
最后编辑时间(Last modified)： 2026/10/21 18:01