R 学习笔记

R 数据类型

向量

概念类似于Python里的列表，只能存储相同类型的数据。

x = 1 # 创建标量
# 创建向量
x <- c(1,2,3)
mode(x)
[1] "numeric"

z <- c(T,F,T)
mode(z)
[1] "logical"

x <- c(1:100) # 1到100的向量
seq(from=1, to=20, by=3) # 创建等差数列，参数可省略。
[1]  1  4  7 10 13 16 19
?seq # 获取帮助

rep(x, 2) # x的每个元素重复2次
[1] 1 2 3 1 2 3
rep(x,c(1,2,3)) # 重复指定次数
[1] 1 2 2 3 3 3

# 向量运算
y <- c(2,3,4)
x+y
[1] 3 5 7
x**2
[1] 1 4 9

向量索引

x <- c(1:100)
length(x) # 输出向量长度
[1] 100

x[1]
[1] 1
x[0] # R从1开始
integer(0)

x[-10] # 不输出第10个
 [1]   1   2   3   4   5   6   7   8   9  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25
[25]  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49
[49]  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73
[73]  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97
[97]  98  99 100

x[-c(10,1)] # 不输出第10个和第1个
 [1]   2   3   4   5   6   7   8   9  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26
[25]  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50
[49]  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74
[73]  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98
[97]  99 100

x[c(2:10)] # 批量索引
[1]  2  3  4  5  6  7  8  9 10
x[c(1, 2, 5, 10)] # 选择索引
[1]  1  2  5 10
x[c(1, 2, 5, 10, 10, 10)] # 可以重复索引
[1]  1  2  5 10 10 10

x[-10, -1]
Error in x[-10, -1] : 量度数目不对 # 矛盾

x[x>3] # 简单的if语句
 [1]   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27
[25]  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
[49]  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75
[73]  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99
[97] 100
x[x>5 & x<8]
[1] 6 7

y <- c(1:10)
y[c(T)] # 逻辑语句
 [1]  1  2  3  4  5  6  7  8  9 10
y[c(F)]
integer(0)
y[c(T,F)] # 两个一循环
[1] 1 3 5 7 9
y[c(T,T)]
 [1]  1  2  3  4  5  6  7  8  9 10
y[c(F,T)]
[1]  2  4  6  8 10
y[c(T,T,T,T,T,T,T,T,T,T,T,T)]
 [1]  1  2  3  4  5  6  7  8  9 10 NA NA # 不足以NA代替

z <- c('a', 'b', 'c', 'd', 'e')
z
[1] "a" "b" "c" "d" "e"
names(z) <- c('one', 'two', 'three', 'four', 'five') # 给z的元素添加名字
z
  one   two three  four  five 
  "a"   "b"   "c"   "d"   "e" 
z['one'] # 索引
one 
"a" 
z['one'] <- 'b' # 赋值
z
  one   two three  four  five 
  "b"   "b"   "c"   "d"   "e" 
z[c(2, 3)] <- c('x', 'y') # 批量赋值
z
  one   two three  four  five 
  "b"   "x"   "y"   "d"   "e" 
z[c(6, 7)] <- c('x', 'y')
z
  one   two three  four  five             
  "b"   "x"   "y"   "d"   "e"   "x"   "y" 
'b' %in% z # 判断
[1] TRUE
z[z %in% c('a', 'b', 'd', 'e')]
 one four five 
 "b" "d"  "e" 
append(x = z, values = 'zzz', after = 0) # 插入元素
        one   two three  four  five             
"zzz"   "b"   "x"   "y"   "d"   "e"   "x"   "y" 
append(x = z, values = 'zzz', after = 1)
  one         two three  four  five             
  "b" "zzz"   "x"   "y"   "d"   "e"   "x"   "y" 
append(x = z, values = 'zzz', after = 2)
  one   two       three  four  five             
  "b"   "x" "zzz"   "y"   "d"   "e"   "x"   "y" 
append(x = z, values = 'zzz', after = 5)
  one   two three  four  five                   
  "b"   "x"   "y"   "d"   "e" "zzz"   "x"   "y" 
append(x = z, values = 'zzz', after = 10)
  one   two three  four  five                   
  "b"   "x"   "y"   "d"   "e"   "x"   "y" "zzz" 
rm(z) # 删除向量
z
错误: 找不到对象'z'

# 删除向量中的元素（不常用）
y
                                      5 
"1" "2" "3" "4" "5" "6" "7" "8" "9" "e" 
y[-c(10)]
                                    
"1" "2" "3" "4" "5" "6" "7" "8" "9" 
y <- y[-c(10)]
y
                                    
"1" "2" "3" "4" "5" "6" "7" "8" "9" 

euro
        ATS         BEF         DEM         ESP         FIM         FRF         IEP         ITL 
  13.760300   40.339900    1.955830  166.386000    5.945730    6.559570    0.787564 1936.270000 
        LUF         NLG         PTE 
  40.339900    2.203710  200.482000 
euro['PTE']
    PTE 
200.482

矩阵与数组

矩阵是二维的，有行和列，每个元素类型要相同，可以看成向量的扩充。

# R内置数据集
head(iris3) # 三种鸢尾花数据
, , Setosa

     Sepal L. Sepal W. Petal L. Petal W.
[1,]      5.1      3.5      1.4      0.2
[2,]      4.9      3.0      1.4      0.2
[3,]      4.7      3.2      1.3      0.2
[4,]      4.6      3.1      1.5      0.2
[5,]      5.0      3.6      1.4      0.2
[6,]      5.4      3.9      1.7      0.4

, , Versicolor

     Sepal L. Sepal W. Petal L. Petal W.
[1,]      7.0      3.2      4.7      1.4
[2,]      6.4      3.2      4.5      1.5
[3,]      6.9      3.1      4.9      1.5
[4,]      5.5      2.3      4.0      1.3
[5,]      6.5      2.8      4.6      1.5
[6,]      5.7      2.8      4.5      1.3

, , Virginica

     Sepal L. Sepal W. Petal L. Petal W.
[1,]      6.3      3.3      6.0      2.5
[2,]      5.8      2.7      5.1      1.9
[3,]      7.1      3.0      5.9      2.1
[4,]      6.3      2.9      5.6      1.8
[5,]      6.5      3.0      5.8      2.2
[6,]      7.6      3.0      6.6      2.1

head(state.x77) # 美国50个州8个指标
           Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
California      21198   5114        1.1    71.71   10.3    62.6    20 156361
Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766

heatmap(state.x77) # 矩阵可直接绘制热图

# 创建矩阵
m <- matrix(1:20, 4,5)
m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
m <- matrix(1:20, 4,5, byrow = T) # 按行分配
m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
m <- matrix(c(1:20), 4,6, byrow = T) # 行列数必须是数据个数的因数
Warning message:
In matrix(c(1:20), 4, 6, byrow = T) : 数据长度[20]不是矩阵列数[6]的整倍数
m <- matrix(c(1:20), 3,5, byrow = T)
Warning message:
In matrix(c(1:20), 3, 5, byrow = T) : 数据长度[20]不是矩阵行数[3]的整倍
m <- matrix(c(1:20), 4,4, byrow = T) # 多时按顺序分配
m
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
[4,]   13   14   15   16
m <- matrix(c(1:20), 4,10, byrow = T) # 少时重复分配
m
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    2    3    4    5    6    7    8    9    10
[2,]   11   12   13   14   15   16   17   18   19    20
[3,]    1    2    3    4    5    6    7    8    9    10
[4,]   11   12   13   14   15   16   17   18   19    20

m <- matrix(1:20, 2,4)
m
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    7
[2,]    2    4    6    8
cnames <- c('C1','C2','C3','C4')
rnames <- c('R1','R2')
dimnames(m) <- list(rnames,cnames) # 命名行列
m
   C1 C2 C3 C4
R1  1  3  5  7
R2  2  4  6  8

dim(m) # 显示对象的维数
[1] 2 4
dim(m) <- c(1,8) # 改变行列数
m
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
dim(m) <- c(4,2)
m
     [,1] [,2]
[1,]    1    5
[2,]    2    6
[3,]    3    7
[4,]    4    8

x <- 1:9
dim(x) <- c(3,3) # 生成矩阵
x
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

矩阵索引

x
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
x[1,1]
[1] 1
x[1,2]
[1] 4
x[1,] # 第1行
[1] 1 4 7
x[,1] # 第1列
[1] 1 2 3
x[1,c(2,3)] # 索引多个
[1] 4 7
x[c(1,2),c(1,2)] # 索引子集
     [,1] [,2]
[1,]    1    4
[2,]    2    5
x[1] # 不加逗号显示行第一个
[1] 1
x[2]
[1] 2
x[-1] #负号表示不显示
[1] 2 3 4 5 6 7 8 9
x[-1,]
     [,1] [,2] [,3]
[1,]    2    5    8
[2,]    3    6    9
x[-1,2]
[1] 5 6
x[1,-2]
[1] 1 7

head(state.x77)
           Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
California      21198   5114        1.1    71.71   10.3    62.6    20 156361
Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766
state.x77['California',] # 通过dimname索引
Population     Income Illiteracy   Life Exp     Murder    HS Grad      Frost       Area 
  21198.00    5114.00       1.10      71.71      10.30      62.60      20.00  156361.00

矩阵运算与线性代数中一致。

m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
n
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    6   10   14   18
[2,]    3    7   11   15   19
[3,]    4    8   12   16   20
[4,]    5    9   13   17   21
m+n
     [,1] [,2] [,3] [,4] [,5]
[1,]    3   11   19   27   35
[2,]    5   13   21   29   37
[3,]    7   15   23   31   39
[4,]    9   17   25   33   41
m*n # 内积，对应元素想乘
     [,1] [,2] [,3] [,4] [,5]
[1,]    2   30   90  182  306
[2,]    6   42  110  210  342
[3,]   12   56  132  240  380
[4,]   20   72  156  272  420
m %*% t(n) # 外积，即线代中的A×B
     [,1] [,2] [,3] [,4]
[1,]  610  655  700  745
[2,]  660  710  760  810
[3,]  710  765  820  875
[4,]  760  820  880  940

rowSums(m) # 行和
[1] 45 50 55 60
colSums(m) # 列和
[1] 10 26 42 58 74
rowMeans(m) # 行平均值
[1]  9 10 11 12
colMeans(m) # 列平均值
[1]  2.5  6.5 10.5 14.5 18.5

n
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
diag(n) # 方阵对角线的值
[1] 1 5 9

数组即多维矩阵。

dim1 <- c('A1', 'A2', 'A3', 'A4', 'A5')
dim2 <- c('B1', 'B2')
dim3 <- c('C1', 'C2')
x <- array(1:20, c(5,2,2), dimnames = list(dim1,dim2,dim3)) # array创建数组
x
, , C1

   B1 B2
A1  1  6
A2  2  7
A3  3  8
A4  4  9
A5  5 10

, , C2

   B1 B2
A1 11 16
A2 12 17
A3 13 18
A4 14 19
A5 15 20

列表

R中最复杂的数据结构，可以用来存储不同的数据，例如向量、矩阵、数据框、列表等。

a <- 1:20
b <- c('a','b','c','d')
mode(mtcars)
[1] "list"
c <- mtcars
d <- matrix(1:20, 4,5)
mlist <- list(a,b,c,d) # 创建列表
mlist
[[1]]
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

[[2]]
[1] "a" "b" "c" "d"

[[3]]
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
...

[[4]]
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

mlist <- list(first=a,second=b,third=c,fourth=d)
mlist
$first
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

$second
[1] "a" "b" "c" "d"

$third
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
...

$fourth
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

mlist <- list('first'=a)
mlist
$first
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

列表索引

mlist['second']
$second
[1] "a" "b" "c" "d"

mlist[2] # 一个[]输出列表的子集
$second
[1] "a" "b" "c" "d"
class(mlist[2])
[1] "list"

mlist[[2]] # 两个[]输出数元素本身的数据类型
[1] "a" "b" "c" "d"
class(mlist[[2]])
[1] "character"
class(mlist$third)
[1] "data.frame"
class(mlist$fourth)
[1] "matrix" "array" 

mlist[5] <- iris # 要用两个[]
Warning message:
In mlist[5] <- iris :
  number of items to replace is not a multiple of replacement length
mlist[[5]] <- iris

# 删除元素
mlist <- mlist[-5]
mlist[4] <- NULL
mlist
$first
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

$second
[1] "a" "b" "c" "d"

$third
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
...

数据框

类似于表格，每一列长度和数据类型必须相同，且必须命名，每行可不相同。

class(iris)
[1] "data.frame"
class(mtcars)
[1] "data.frame"
iris
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
...
mtcars
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
...

# 创建数据框
state <- data.frame(state.name,state.abb,state.x77)
state
                   state.name state.abb Population Income Illiteracy Life.Exp Murder HS.Grad Frost
Alabama               Alabama        AL       3615   3624        2.1    69.05   15.1    41.3    20
Alaska                 Alaska        AK        365   6315        1.5    69.31   11.3    66.7   152
Arizona               Arizona        AZ       2212   4530        1.8    70.55    7.8    58.1    15
...
rownames(state)
 [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"       "California"    
 [6] "Colorado"       "Connecticut"    "Delaware"       "Florida"        "Georgia"       
[11] "Hawaii"         "Idaho"          "Illinois"       "Indiana"        "Iowa"          
[16] "Kansas"         "Kentucky"       "Louisiana"      "Maine"          "Maryland"      
[21] "Massachusetts"  "Michigan"       "Minnesota"      "Mississippi"    "Missouri"      
[26] "Montana"        "Nebraska"       "Nevada"         "New Hampshire"  "New Jersey"    
[31] "New Mexico"     "New York"       "North Carolina" "North Dakota"   "Ohio"          
[36] "Oklahoma"       "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina"
[41] "South Dakota"   "Tennessee"      "Texas"          "Utah"           "Vermont"       
[46] "Virginia"       "Washington"     "West Virginia"  "Wisconsin"      "Wyoming"       
colnames(state)
 [1] "state.name" "state.abb"  "Population" "Income"     "Illiteracy" "Life.Exp"   "Murder"    
 [8] "HS.Grad"    "Frost"      "Area" 

# 索引
state[1]
                   state.name
Alabama               Alabama
Alaska                 Alaska
Arizona               Arizona
...
state['state.abb']
               state.abb
Alabama               AL
Alaska                AK
Arizona               AZ
...
state$Murder
 [1] 15.1 11.3  7.8 10.1 10.3  6.8  3.1  6.2 10.7 13.9  6.2  5.3 10.3  7.1  2.3  4.5 10.6 13.2  2.7
[20]  8.5  3.3 11.1  2.3 12.5  9.3  5.0  2.9 11.5  3.3  5.2  9.7 10.9 11.1  1.4  7.4  6.4  4.2  6.1
[39]  2.4 11.6  1.7 11.0 12.2  4.5  5.5  9.5  4.3  6.7  3.0  6.9
state[,'state.abb']
 [1] "AL" "AK" "AZ" "AR" "CA" "CO" "CT" "DE" "FL" "GA" "HI" "ID" "IL" "IN" "IA" "KS" "KY" "LA" "ME"
[20] "MD" "MA" "MI" "MN" "MS" "MO" "MT" "NE" "NV" "NH" "NJ" "NM" "NY" "NC" "ND" "OH" "OK" "OR" "PA"
[39] "RI" "SC" "SD" "TN" "TX" "UT" "VT" "VA" "WA" "WV" "WI" "WY"
state['California',]
           state.name state.abb Population Income Illiteracy Life.Exp Murder HS.Grad Frost   Area
California California        CA      21198   5114        1.1    71.71   10.3    62.6    20 156361

women
   height weight
1      58    115
2      59    117
3      60    120
...
class(state['state.abb']) # 单括号仍为数据框
[1] "data.frame"
class(state[['state.abb']]) # 双括号为数据本身
[1] "character"
state[['state.abb']]
 [1] "AL" "AK" "AZ" "AR" "CA" "CO" "CT" "DE" "FL" "GA" "HI" "ID" "IL" "IN" "IA" "KS" "KY" "LA" "ME"
[20] "MD" "MA" "MI" "MN" "MS" "MO" "MT" "NE" "NV" "NH" "NJ" "NM" "NY" "NC" "ND" "OH" "OK" "OR" "PA"
[39] "RI" "SC" "SD" "TN" "TX" "UT" "VT" "VA" "WA" "WV" "WI" "WY"
state['state.abb']
               state.abb
Alabama               AL
Alaska                AK
Arizona               AZ
...

因子

因子是R中重要的概念。

R中变量类型

名义型变量：没有顺序分别，相互独立，如省份名。
有序型变量：介于二者之间，不同值之间有顺序关系，但不是连续的数量变化，如good、better、best。
连续型变量：某个范围中的任意值，如存储年龄、身高、GDP的变量。

一般，字符型数据多为名义型变量，数值型数据多为连续型变量。

名义型变量和有序型变量在R中称为因子，这些变量的可能值称为一个水平level，如good，better，best都称为一个level。

因子的最大作用是用来分类，计算频数和频率。

# 内置数据集
state.region
 [1] South         West          West          South         West          West         
 [7] Northeast     South         South         South         West          West         
[13] North Central North Central North Central North Central South         South        
...         
Levels: Northeast South North Central West

state.division
 [1] East South Central Pacific            Mountain           West South Central
 [5] Pacific            Mountain           New England        South Atlantic    
 [9] South Atlantic     South Atlantic     Pacific            Mountain          
...        
9 Levels: New England Middle Atlantic South Atlantic ... Pacific

# 定义因子
f <- factor(c(1,1,2,2,2,3,4))
f
[1] 1 1 2 2 2 3 4
Levels: 1 2 3 4

mtcars$cyl
 [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4
factor(mtcars$cyl)
 [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4
Levels: 4 6 8

# 可以设置顺序
factor(c('one','one', 'two', 'three', 'three'), ordered = T, levels = c('one', 'two', 'three', 'four'))
[1] one   one   two   three three
Levels: one < two < three < four
factor(c('five','five', 'two', 'three', 'three'), ordered = T, levels = c('one', 'two', 'three', 'four'))
[1] <NA <NA two   three three # 因子不能含有没有的level
Levels: one < two < three < four

plot(mtcars$cyl) # 散点图
plot(factor(mtcars$cyl)) # 因子做出条形图

# cut函数可将向量按要求分隔为因子
n <- 1:20
cut(n,c(seq(0,20,5)))
 [1] (0,5]   (0,5]   (0,5]   (0,5]   (0,5]   (5,10]  (5,10]  (5,10]  (5,10]  (5,10]  (10,15]
[12] (10,15] (10,15] (10,15] (10,15] (15,20] (15,20] (15,20] (15,20] (15,20]
Levels: (0,5] (5,10] (10,15] (15,20]

mtcars$wt
 [1] 2.620 2.875 2.320 3.215 3.440 3.460 3.570 3.190 3.150 3.440 3.440 4.070 3.730 3.780 5.250
[16] 5.424 5.345 2.200 1.615 1.835 2.465 3.520 3.435 3.840 3.845 1.935 2.140 1.513 3.170 2.770
[31] 3.570 2.780
range(mtcars$wt)
[1] 1.513 5.424
plot(cut(mtcars$wt,c(seq(1,6,1)))) # 直接画出频数条形图

缺失数据

NA（not available）代表缺失值，用来存储缺失数据，表示没有。注意没有不是0，可能是任何值。

NaN代表不可能值，不存在。

Inf表示无穷大，-Inf为无穷小。

在R中有相关包和函数可以处理缺失值，如VIM包。

[1] NaN
NA+1
[1] NA
0/0
[1] NaN
1/0
[1] Inf
-1/0
[1] -Inf

a <- c(NA, NaN, Inf, -Inf)
is.na(a)
[1]  TRUE  TRUE FALSE FALSE
is.infinite(a)
[1] FALSE FALSE  TRUE  TRUE
is.nan(a)
[1] FALSE  TRUE FALSE FALSE

a
 [1] NA  1  2  3  4  5  6  7  8  9 10
sum(a)
[1] NA
mean(a)
[1] NA
sum(a, na.rm = T) # na.rm = T，计算时跳过缺失值
[1] 55
mean(a, na.rm = T)
[1] 5.5

a
[1] NA  1  2  3  4 NA NA
b <- na.omit(a) # 去除缺失数据
b
[1] 1 2 3 4
attr(,"na.action")
[1] 1 6 7
attr(,"class")
[1] "omit"
is.na(b)
[1] FALSE FALSE FALSE FALSE
sum(b)
[1] 10

a
       state.name state.area
1            <NA     51609
2            <NA    589757
3         Arizona     113909
4        Arkansas      53104
5      California     158693

na.omit(a) # 为数据框时删除此行
       state.name state.area
3         Arizona     113909
4        Arkansas      53104
5      California     158693

字符串

# 计算字符串长度
nchar('hello')
[1] 5
a <- c('hell o', 'world', '!') # 空格也算
nchar(a)
[1] 6 5 1
month.name
 [1] "January"   "February"  "March"     "April"     "May"       "June"      "July"     
 [8] "August"    "September" "October"   "November"  "December" 
nchar(month.name)
 [1] 7 8 5 5 3 4 4 6 9 7 8 8
length(month.name) # 向量的长度，注意区分二者
[1] 12
nchar(12345) # 为数字时当成字符串处理
[1] 5
nchar(c(12, 133))
[1] 2 3

# 字符串拼接
paste('hello', 'world', '!')
[1] "hello world !"
paste('hello', 'world', '!', sep = '-') # 可选分隔符，默认空格
[1] "hello-world-!"
a <- c('Zhao', 'Qian', 'Sun') # 注意向量
paste(a, 'Ming')
[1] "Zhao Ming" "Qian Ming" "Sun Ming" 

# 字符串提取
substr(x = month.name, start = 1, stop = 3)
 [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"

# 大小写转换
toupper(month.name)
 [1] "JANUARY"   "FEBRUARY"  "MARCH"     "APRIL"     "MAY"       "JUNE"      "JULY"     
 [8] "AUGUST"    "SEPTEMBER" "OCTOBER"   "NOVEMBER"  "DECEMBER" 
tolower(month.name)
 [1] "january"   "february"  "march"     "april"     "may"       "june"      "july"     
 [8] "august"    "september" "october"   "november"  "december" 

# 局部大小写转换-正则表达式
# 待学
gsub()

# 字符串匹配
a <- c('Zhao', 'Qian', 'Sun', 'Li')
grep('Sun', a)
[1] 3
grep('a', a)
[1] 1 2
grep('ao', a)
[1] 1
grep('au', a)
integer(0)
a <- c('A+', 'AC')
grep('A+', a, fixed = T)
[1] 1
grep('A+', a, fixed = F) # 含有fixed参数表示支持正则表达式，+表示任何字符
[1] 1 2
match('A', a) # match不支持正则表达式
[1] NA
match('A+', a)
[1] 1
match('AC', a)
[1] 2

# 字符串分隔
strsplit(x = 'hello my dear.', split = ' ')
[[1]]
[1] "hello" "my"    "dear."

a <- c('oh my god!', 'hello my dear.')
strsplit(a, ' ') # 可分隔向量，返回list
[[1]]
[1] "oh"   "my"   "god!"

[[2]]
[1] "hello" "my"    "dear."

# 生成扑克牌
a
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13
b
[1] "heitao"   "fangpian" "hongtao"  "meihua"  
outer(b, a, FUN = paste)
     [,1]         [,2]         [,3]         [,4]         [,5]         [,6]        
[1,] "heitao 1"   "heitao 2"   "heitao 3"   "heitao 4"   "heitao 5"   "heitao 6"  
[2,] "fangpian 1" "fangpian 2" "fangpian 3" "fangpian 4" "fangpian 5" "fangpian 6"
[3,] "hongtao 1"  "hongtao 2"  "hongtao 3"  "hongtao 4"  "hongtao 5"  "hongtao 6" 
[4,] "meihua 1"   "meihua 2"   "meihua 3"   "meihua 4"   "meihua 5"   "meihua 6"  
     [,7]         [,8]         [,9]         [,10]         [,11]         [,12]        
[1,] "heitao 7"   "heitao 8"   "heitao 9"   "heitao 10"   "heitao 11"   "heitao 12"  
[2,] "fangpian 7" "fangpian 8" "fangpian 9" "fangpian 10" "fangpian 11" "fangpian 12"
[3,] "hongtao 7"  "hongtao 8"  "hongtao 9"  "hongtao 10"  "hongtao 11"  "hongtao 12" 
[4,] "meihua 7"   "meihua 8"   "meihua 9"   "meihua 10"   "meihua 11"   "meihua 12"  
     [,13]        
[1,] "heitao 13"  
[2,] "fangpian 13"
[3,] "hongtao 13" 
[4,] "meihua 13"
outer(b, a, FUN = paste, sep = '-')

正则表达式

regular expression regex RE

特点

正则表达式能够凝练字符串的特征
是通用的字符串表达框架
是简洁表达一组字符串的表达式
可以判断某个字符串的特征归属，是否具有这个特征

作用

表达文本类型的特征

同时查找或替换一组字符串

匹配字符串的全部或部分

语法

由字符和操作符构成

基本和常用的正则表达式操作符如下

操作符	说明	例子
`.`	表示任何单个字符
`[]`	给出单个字符的取值范围	`[abc]`：a或b或c，`[a-z]`：a-z单个字母
`[^]`	排除此取值范围	`[^abc]`：非a或b或c
`*`	前一个字符出现`>=0`次	`abc*`：ab、abc、abcc、abccc等
`+`	前一个字符出现`>=1`次	`abc+`：abc、abcc、abccc等
`?`	前一个字符出现`0`或`1`次	`abc?`：ab、abc
`	`	左右表达式中的一个
`{m}`	前一个字符出现`m`次	`ab{2}c`：abbc
`{n,}`	前一个字符出现`>=n`次	`ab{2,}c`：abbc、abbbc、abbbbc等
`{m,n}`	前一个字符出现`m-n`次（含`n`）	`ab{1,3}c`：abc、abbc、abbbc
`^`	匹配字符串的开头	`^123`：字符串以abc开头
`$`	匹配字符串的结尾	`123$`：字符串以123结尾
`()`	分组标记，内部只能使用`	`
`\d`	数字，等价于`[0-9]`
`\w`	单词字符，等价于`[A-Za-z0-9_]`

经典实例

^[A-Za-z]+$：由字母组成的字符串

^[A-Za-z]+$：由字母和数字组成的字符串

^-?[1-9]\d*$：整数

^[1-9]\d*$：正整数

[\u4e00-\u9fa5]：匹配中文字符

日期和时间

时间序列数据

# 内置时间序列数据
class(sunspots)
[1] "ts"
sunspots
       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
1749  58.0  62.6  70.0  55.7  85.0  83.5  94.8  66.3  75.9  75.5 158.6  85.2
1750  73.3  75.9  89.2  88.3  90.0 100.0  85.4 103.0  91.2  65.7  63.3  75.4
1751  70.0  43.5  45.3  56.4  60.7  50.7  66.3  59.8  23.5  23.2  28.5  44.0
...

presidents
     Qtr1 Qtr2 Qtr3 Qtr4
1945   NA   87   82   75
1946   63   50   43   32
1947   35   60   54   55
...
airmiles
Time Series:
Start = 1937 
End = 1960 
Frequency = 1 
 [1]   412   480   683  1052  1385  1418  1634  2178  3362  5948  6109  5981  6753  8003 10566
[16] 12528 14760 16769 19819 22362 25340 25343 29269 30514
class(airmiles)
[1] "ts"

Sys.Date() # 输出系统时间
[1] "2022-01-31"
class(Sys.Date())
[1] "Date" # 单独的数据类型，不是字符串

a <- as.Date('2022-02-02', format = '%Y-%m-%d') # 字符串转换成时间
class(a)
[1] "Date"
?strftime # 更多格式信息查看

# 生成时间序列，系统自动按时间处理，不会出现2月30日等情况
seq(as.Date('2022-01-01'), as.Date('2022-03-01'), by=5)
 [1] "2022-01-01" "2022-01-06" "2022-01-11" "2022-01-16" "2022-01-21" "2022-01-26" "2022-01-31"
 [8] "2022-02-05" "2022-02-10" "2022-02-15" "2022-02-20" "2022-02-25"

# 生成随机数
runif(10, min = 0, max = 10)
 [1] 9.0326531 0.5515827 6.8552402 0.8496375 7.1430260 5.6753642 9.7609197 7.2243469 5.2854405
[10] 7.4681394

sales <- round(runif(48, min = 50, max = 100), 0)
sales
 [1] 98 96 77 51 55 54 80 71 59 65 70 81 54 75 64 99 74 58 89 93 78 80 83 71 53 99 69 64 91 94
[31] 72 54 85 83 58 74 55 76 54 57 87 79 83 72 54 87 96 61

ts(sales, start = c(2022,1), end = c(2022, 12), frequency = 12) # 生成时间序列
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2022  98  96  77  51  55  54  80  71  59  65  70  81

# frequency：1表示按年，12表示按月份，4表示按季度
ts(sales, start = c(2022,1), end = c(2025, 1), frequency = 1)
Time Series:
Start = 2022 
End = 2025 
Frequency = 1 
[1] 98 96 77 51
ts(sales, start = c(2022,1), end = c(2023, 1), frequency = 12)
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2022  98  96  77  51  55  54  80  71  59  65  70  81
2023  54                                            
ts(sales, start = c(2022,1), end = c(2023, 1), frequency = 4)
     Qtr1 Qtr2 Qtr3 Qtr4
2022   98   96   77   51
2023   55

数据处理

数据获取

获取数据的途径：

键盘直接敲入
读取存储在外部文件中已预处理好的数据
访问数据库系统

# 通过键盘敲入
patientID <- c(1, 2, 3, 4)
admdate <- c("10/15/2009", "11/01/2009", "10/21/2009", "10/28/2009")
age <- c(25, 34, 28, 52)
diabetes <- c("Type1", "Type2", "Type1", "Type1")
status <- c("Poor", "Improved", "Excellent", "Poor")
data <- data.frame(patientID, age, diabetes, status)

# 或通过编辑器输入，需先设置变量，否则数据无法保存
data2 <- data.frame(patientID=character(), admdate=character(), age=numeric(), diabetes=character(), status=character())
data2 <- edit(data2)
# 或fix()函数修改，不用再赋值
fix(data2)

# 访问数据库系统
install.packages("RODBC")

文件读取

纯文本文件.csv .txt

read.table('/home/lah/Desktop/R_DATA/input.txt')
    Ozone Solar.R Wind Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
...

# 默认空格分隔符
x <- read.table('/home/lah/Desktop/R_DATA/input.csv')
x
                    V1                                                                  V2
1                      ,"mpg","cyl","disp","hp","drat","wt","qsec","vs","am","gear","carb"
2            Mazda RX4                                ,21,6,160,110,3.9,2.62,16.46,0,1,4,4
3        Mazda RX4 Wag                               ,21,6,160,110,3.9,2.875,17.02,0,1,4,4
...

x <- read.table('/home/lah/Desktop/R_DATA/input.csv', sep = ',')
head(x)
                 V1   V2  V3   V4  V5   V6    V7    V8 V9 V10  V11  V12
1                    mpg cyl disp  hp drat    wt  qsec vs  am gear carb
2         Mazda RX4   21   6  160 110  3.9  2.62 16.46  0   1    4    4
3     Mazda RX4 Wag   21   6  160 110  3.9 2.875 17.02  0   1    4    4
...

x <- read.table('/home/lah/Desktop/R_DATA/input.csv', sep = ',', header = T, row.names = 1)
head(x)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

# 可跳过前多少行
x <- read.table('/home/lah/Desktop/R_DATA/input.csv', sep = ',', skip = 30)
x
             V1   V2 V3  V4  V5   V6   V7   V8 V9 V10 V11 V12
1  Ferrari Dino 19.7  6 145 175 3.62 2.77 15.5  0   1   5   6
2 Maserati Bora 15.0  8 301 335 3.54 3.57 14.6  0   1   5   8
3    Volvo 142E 21.4  4 121 109 4.11 2.78 18.6  1   1   4   2

# 只读取多少行
x <- read.table('/home/lah/Desktop/R_DATA/input.csv', sep = ',', nrows = 2)
x
         V1  V2  V3   V4  V5   V6   V7    V8 V9 V10  V11  V12
1           mpg cyl disp  hp drat   wt  qsec vs  am gear carb
2 Mazda RX4  21   6  160 110  3.9 2.62 16.46  0   1    4    4

read.fwf('/home/lah/Desktop/R_DATA/fwf.txt',widths = c(11,11))

# read.table()的简化版本
read.csv(file, header = TRUE, sep = ",", quote = "\"",
         dec = ".", fill = TRUE, comment.char = "", ...)

read.csv2(file, header = TRUE, sep = ";", quote = "\"",
          dec = ",", fill = TRUE, comment.char = "", ...)

read.delim(file, header = TRUE, sep = "\t", quote = "\"",
           dec = ".", fill = TRUE, comment.char = "", ...)

read.delim2(file, header = TRUE, sep = "\t", quote = "\"",
            dec = ",", fill = TRUE, comment.char = "", ...)

网页表格，剪贴板，压缩文件

x <- read.table("https://codeload.github.com/mperdeck/LINQtoCSV/zip/master", header = TRUE)

library(XML)
readHTMLTable("https://en.wikipedia.org/wiki/World_population", which=3)

# 搜索相关R包
RSiteSearch('manhattan')
A search query has been submitted to http://search.r-project.org
计算结果应很快就在瀏覽器里打开

# 读取剪贴板
read.table('clipboard', header = T)

# 直接读取压缩文件，不用解压
read.table(gzfile('/home/lah/Desktop/R_DATA/newfile.txt.gz'))
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
...

scan()和readLines()

readLines('/home/lah/Desktop/R_DATA/newfile.txt', n=4)
[1] "\"Ozone\" \"Solar.R\" \"Wind\" \"Temp\" \"Month\" \"Day\""
[2] "\"1\" 41 190 7.4 67 5 1"                                  
[3] "\"2\" 36 118 8 72 5 2"                                    
[4] "\"3\" 12 149 12.6 74 5 3" 

readLines('/home/lah/Desktop/R_DATA/scan.txt')
 [1] "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6"
 [4] "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6"
 [7] "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6" "one\t2\t3\tfour\t5\t6"
[10] "one\t2\t3\tfour\t5\t6" "" 

# scan()按特定方式扫描文件
read.table('/home/lah/Desktop/R_DATA/scan.txt', sep = '\t')
    V1 V2 V3   V4 V5 V6
1  one  2  3 four  5  6
2  one  2  3 four  5  6
3  one  2  3 four  5  6
4  one  2  3 four  5  6
5  one  2  3 four  5  6
6  one  2  3 four  5  6
7  one  2  3 four  5  6
8  one  2  3 four  5  6
9  one  2  3 four  5  6
10 one  2  3 four  5  6

scan('/home/lah/Desktop/R_DATA/scan.txt', what = list(character(3)))
Read 60 records
[[1]]
 [1] "one"  "2"    "3"    "four" "5"    "6"    "one"  "2"    "3"    "four" "5"    "6"    "one" 
[14] "2"    "3"    "four" "5"    "6"    "one"  "2"    "3"    "four" "5"    "6"    "one"  "2"   
[27] "3"    "four" "5"    "6"    "one"  "2"    "3"    "four" "5"    "6"    "one"  "2"    "3"   
[40] "four" "5"    "6"    "one"  "2"    "3"    "four" "5"    "6"    "one"  "2"    "3"    "four"
[53] "5"    "6"    "one"  "2"    "3"    "four" "5"    "6" 

scan('/home/lah/Desktop/R_DATA/scan.txt', what = list(character(3), numeric(1), numeric(1)))
Read 20 records
[[1]]
 [1] "one"  "four" "one"  "four" "one"  "four" "one"  "four" "one"  "four" "one"  "four" "one" 
[14] "four" "one"  "four" "one"  "four" "one"  "four"

[[2]]
 [1] 2 5 2 5 2 5 2 5 2 5 2 5 2 5 2 5 2 5 2 5

[[3]]
 [1] 3 6 3 6 3 6 3 6 3 6 3 6 3 6 3 6 3 6 3 6

scan('/home/lah/Desktop/R_DATA/scan.txt', what = list(A=character(3), B=numeric(1), C=numeric(1), D=character(3), E=numeric(1), F=numeric(1)))
Read 10 records
$A
 [1] "one" "one" "one" "one" "one" "one" "one" "one" "one" "one"

$B
 [1] 2 2 2 2 2 2 2 2 2 2

$C
 [1] 3 3 3 3 3 3 3 3 3 3

$D
 [1] "four" "four" "four" "four" "four" "four" "four" "four" "four" "four"

$E
 [1] 5 5 5 5 5 5 5 5 5 5

$F
 [1] 6 6 6 6 6 6 6 6 6 6

文件写入

# write()直接写入向量或矩阵
write(x, file = '2.txt')

# 一行写入的个数
write(x, file = '2.txt', ncolumns = 10)
write(x, file = '2.txt', ncolumns = 10, sep = '\t')

# write.table()写入数据框或矩阵
write.table(state.x77, '2.txt')
write.table(mtcars, '2.txt')

# 设置分隔符
write.table(mtcars, '2.csv', sep = ',')

# 追加
write.table(mtcars, '2.csv', append = T)

# 字符串不加引号
write.table(mtcars, '2.csv', quote = F)

# 去除行名或列名
write.table(mtcars, '2.csv', row.names = F)
write.table(mtcars, '2.csv', col.names = F)

# 直接写成压缩文件
write.table(mtcars, gzfile('1.txt.gz'))

# 写入或读取其他统计软件格式的文件，用foreign包
library(foreign)
help(package = 'foreign')

Excel文件读写

XLConnect包（需Java：https://blog.csdn.net/weixin_36350581/article/details/114214756）运行失败

XLSX包运行失败

openxlsx包

# 读取，默认读取第一个sheet
read.xlsx('writeXLSX2.xlsx')
read.xlsx('writeXLSX2.xlsx', sheet = 'state info')

# 从某行开始读
read.xlsx('writeXLSX2.xlsx', startRow = 48)

# 是否包括行名和列名
read.xlsx('writeXLSX2.xlsx', colNames = F)
read.xlsx('writeXLSX2.xlsx', rowNames = F)

# 是否跳过空行和空列
read.xlsx('writeXLSX2.xlsx', skipEmptyRows = F)
read.xlsx('writeXLSX2.xlsx', skipEmptyCols = F)

# 读取指定行列数
read.xlsx('writeXLSX2.xlsx', sheet = 'test', cols = c(1))
  A
1 1
2 2
3 3
4 4
5 5
6 6
read.xlsx('writeXLSX2.xlsx', sheet = 'test', cols = c(1:4))
  A B C
1 1 2 3
2 2 3 4
3 3 4 5
4 4 5 6
5 5 6 7
6 6 7 8
read.xlsx('writeXLSX2.xlsx', sheet = 'test', cols = c(1:4), rows = c(1:4))
  A B C
1 1 2 3
2 2 3 4
3 3 4 5

# 写入
write.xlsx(mtcars, 'test.xlsx')

# 作为表格写入
write.xlsx(mtcars, 'test.xlsx', asTable = T)

# 默认覆盖写入
write.xlsx(mtcars, 'test.xlsx', overwrite = F)
Error in saveWorkbook(wb, file = file, overwrite = overwrite) : 
  File already exists!

# 写入行名或列名
write.xlsx(mtcars, file = 'car.xlsx', rowNames=T, colNames=T)

# 可写入多个data frame到不同的sheet
l <- list('Cars Info' = mtcars, 'Iris Info' = iris, 'State info' = data.frame(state.name, state.abb, state.x77))
write.xlsx(l, 'test.xlsx')

R 格式文件读写

.rds文件可以保存单个R对象

.rdata文件可以保存对个R对象

存储为R文件有很多优势，特别是对于大数据文件，R会自动压缩处理并包含更多信息。

head(iris, n = 3)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa

# 保存为rds文件
saveRDS(iris, file='iris.rds')

# 读取rds文件为一个R对象，交互界面下双击文件即可
iris_test <- readRDS("~/R_learn/R_Data_Structure/iris.rds")
iris_test
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
...

# 将R对象保存至Rdata文件
save(mtcars, iris, iris3, file = 'test_data.Rdata')

# 载入Rdata文件
# 注意如果要载入的R对象名称与现有工作空间内R对象重名，将覆盖
load("~/R_learn/R_Data_Structure/test_data.Rdata")

# 保存所有R对象
save.image(file = 'test.Rdata')

数据转换

数据类型转换

head(car， n = 3)
                     mpg cyl disp  hp drat    wt  qsec vs am gear carb
1         Mazda RX4 21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
2     Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
3        Datsun 710 22.8   4  108  93 3.85 2.320 18.61  1  1    4    1

# is多在函数编程时作为条件判断
is.data.frame(car)
[1] TRUE
is.matrix(car)
[1] FALSE
is.na(car)
                      mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
Mazda RX4           FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
Mazda RX4 Wag       FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
Datsun 710          FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...

# 查看所有is
methods(is)
 [1] is.array                is.atomic               is.call                
 [4] is.character            is.complex              is.data.frame          
 [7] is.double               is.element              is.empty.model  
...

# as转换数类型
# as并不总是奏效，如数据框无法转换成向量
dstate.x77 <- as.data.frame(state.x77)
View(dstate.x77)
is.data.frame(dstate.x77)
[1] TRUE

# 查看所有as
methods(as)
  [1] as.array                      as.array.default              as.call                      
  [4] as.character                  as.character.condition        as.character.Date            
  [7] as.character.default          as.character.error            as.character.factor 
...

# 向量可以很方便地做各种转换
a <- state.abb
head(a)
[1] "AL" "AK" "AZ" "AR" "CA" "CO"

# 加上维度转换成矩阵
dim(a) <- c(5,10)
head(a, n = 3)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] "AL" "CO" "HI" "KS" "MA" "MT" "NM" "OK" "SD" "VA" 
[2,] "AK" "CT" "ID" "KY" "MI" "NE" "NY" "OR" "TN" "WA" 
[3,] "AZ" "DE" "IL" "LA" "MN" "NV" "NC" "PA" "TX" "WV" 
a <- state.abb

# 转换为因子
as.factor(a)
 [1] AL AK AZ AR CA CO CT DE FL GA HI ID IL IN IA KS KY LA ME MD MA MI MN MS MO MT NE NV NH NJ
[31] NM NY NC ND OH OK OR PA RI SC SD TN TX UT VT VA WA WV WI WY
50 Levels: AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS MT ... WY
as.list(a)
[[1]]
[1] "AL"

[[2]]
[1] "AK"

[[3]]
[1] "AZ"
...

# 组成数据框
state <- data.frame(a,state.x77)
head(state, n = 3)
            a Population Income Illiteracy Life.Exp Murder HS.Grad Frost   Area
Alabama    AL       3615   3624        2.1    69.05   15.1    41.3    20  50708
Alaska     AK        365   6315        1.5    69.31   11.3    66.7   152 566432
Arizona    AZ       2212   4530        1.8    70.55    7.8    58.1    15 113417
state$Population
 [1]  3615   365  2212  2110 21198  2541  3100   579  8277  4931   868   813 11197  5313  2861
[16]  2280  3387  3806  1058  4122  5814  9111  3921  2341  4767   746  1544   590   812  7333
[31]  1144 18076  5441   637 10735  2715  2284 11860   931  2816   681  4173 12237  1203   472
[46]  4981  3559  1799  4589   376

b <- state['California',]
is.data.frame(b)
[1] TRUE

# 不显示名字
unname(b)
                                                      
California CA 21198 5114 1.1 71.71 10.3 62.6 20 156361
head(b)
            a Population Income Illiteracy Life.Exp Murder HS.Grad Frost   Area
California CA      21198   5114        1.1    71.71   10.3    62.6    20 156361

# 数据框转换为向量
unlist(b)
         a Population     Income Illiteracy   Life.Exp     Murder    HS.Grad      Frost 
      "CA"    "21198"     "5114"      "1.1"    "71.71"     "10.3"     "62.6"       "20" 
      Area 
  "156361" 

# 这样不能转换
as.vector(b)
            a Population Income Illiteracy Life.Exp Murder HS.Grad Frost   Area
California CA      21198   5114        1.1    71.71   10.3    62.6    20 156361
c <- as.vector(b)
class(c)
[1] "data.frame"

数据取子集、合并、添加

# 索引
state.x77[1:5, 1:5]
           Population Income Illiteracy Life Exp Murder
Alabama          3615   3624        2.1    69.05   15.1
Alaska            365   6315        1.5    69.31   11.3
Arizona          2212   4530        1.8    70.55    7.8
Arkansas         2110   3378        1.9    70.66   10.1
California      21198   5114        1.1    71.71   10.3

state.x77[c(1,2,3), 1:2]
        Population Income
Alabama       3615   3624
Alaska         365   6315
Arizona       2212   4530

state.x77[c(1,2,7), 1:2]
            Population Income
Alabama           3615   3624
Alaska             365   6315
Connecticut       3100   5348

state.x77[1:3, -1:-5]
        HS Grad Frost   Area
Alabama    41.3    20  50708
Alaska     66.7   152 566432
Arizona    58.1    15 113417

# 条件索引
mtcars[mtcars$cyl==4,]
                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Datsun 710     22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230       22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
...

mtcars[mtcars$mpg>30,]
                mpg cyl disp  hp drat    wt  qsec vs am gear carb
Fiat 128       32.4   4 78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic    30.4   4 75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla 33.9   4 71.1  65 4.22 1.835 19.90  1  1    4    1
Lotus Europa   30.4   4 95.1 113 3.77 1.513 16.90  1  1    5    2

subset(mtcars, mtcars$mpg>32)
                mpg cyl disp hp drat    wt  qsec vs am gear carb
Fiat 128       32.4   4 78.7 66 4.08 2.200 19.47  1  1    4    1
Toyota Corolla 33.9   4 71.1 65 4.22 1.835 19.90  1  1    4    1

mtcars[mtcars$mpg>30 & mtcars$wt>2,]
          mpg cyl disp hp drat  wt  qsec vs am gear carb
Fiat 128 32.4   4 78.7 66 4.08 2.2 19.47  1  1    4    1

# 使用subset函数
subset(mtcars, mtcars$mpg>32 & mtcars$wt>2)
          mpg cyl disp hp drat  wt  qsec vs am gear carb
Fiat 128 32.4   4 78.7 66 4.08 2.2 19.47  1  1    4    1

# sample函数抽样
x <- 1:100
sample(x, 10)
 [1] 46 29 39 20 12 70 66 74 33 94

# 有放回抽样，可以重复
sample(x, 10, replace = T)
 [1] 23 23 22 52 37 20 78 96 34 92

# 对数据框进行抽样
sample(mtcars$mpg, 5)
[1] 22.8 15.2 10.4 21.0 21.0
mtcars[sample(mtcars$mpg, 5), ]
                    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Honda Civic        30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Chrysler Imperial  14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
Toyota Corona      21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Cadillac Fleetwood 10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Fiat X1-9          27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1

# 删除数据框的某些行和列
mtcars$mpg <- NULL
head(mtcars)
                  cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant             6  225 105 2.76 3.460 20.22  1  0    3    1

mtcars[-1:-30,]
              cyl disp  hp drat   wt qsec vs am gear carb
Maserati Bora   8  301 335 3.54 3.57 14.6  0  1    5    8
Volvo 142E      4  121 109 4.11 2.78 18.6  1  1    4    2

mtcars[,-1:-8]
                    gear carb
Mazda RX4              4    4
Mazda RX4 Wag          4    4
Datsun 710             4    1
...

# 数据框的添加与合并
state.abb
 [1] "AL" "AK" "AZ" "AR" "CA" "CO" "CT" "DE" "FL" "GA" "HI" "ID" "IL" "IN" "IA" "KS" "KY" "LA" "ME"
[20] "MD" "MA" "MI" "MN" "MS" "MO" "MT" "NE" "NV" "NH" "NJ" "NM" "NY" "NC" "ND" "OH" "OK" "OR" "PA"
[39] "RI" "SC" "SD" "TN" "TX" "UT" "VT" "VA" "WA" "WV" "WI" "WY"

# 合并列
data.frame(USArrests, state.abb)
               Murder Assault UrbanPop Rape state.abb
Alabama          13.2     236       58 21.2        AL
Alaska           10.0     263       48 44.5        AK
Arizona           8.1     294       80 31.0        AZ
...
   
cbind(USArrests, state.abb)
               Murder Assault UrbanPop Rape state.abb
Alabama          13.2     236       58 21.2        AL
Alaska           10.0     263       48 44.5        AK
Arizona           8.1     294       80 31.0        AZ
...

# 合并行
data1 <- head(mtcars, 5)
data2 <- tail(mtcars, 5)
rbind(data1, data2)
                  cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Lotus Europa        4  95.1 113 3.77 1.513 16.90  1  1    5    2
Ford Pantera L      8 351.0 264 4.22 3.170 14.50  0  1    5    4
Ferrari Dino        6 145.0 175 3.62 2.770 15.50  0  1    5    6
Maserati Bora       8 301.0 335 3.54 3.570 14.60  0  1    5    8
Volvo 142E          4 121.0 109 4.11 2.780 18.60  1  1    4    2

# 重复数据不会自动去除
data2 <- head(mtcars, 10)
rbind(data1, data2)
                   cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4            6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag        6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710           4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive       6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout    8 360.0 175 3.15 3.440 17.02  0  0    3    2
Mazda RX41           6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag1       6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 7101          4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive1      6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout1   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant              6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360           8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D            4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230             4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280             6 167.6 123 3.92 3.440 18.30  1  0    4    4

data3 <- rbind(data1, data2)
# duplicated函数判断是否存在重复数据
duplicated(data3)
 [1] FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE

# 取出重复数据
data3[duplicated(data3), ]
                   cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX41           6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag1       6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 7101          4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive1      6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout1   8  360 175 3.15 3.440 17.02  0  0    3    2

#!对TRUE和FALSE取反
data3[!duplicated(data3), ]
                  cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            6 167.6 123 3.92 3.440 18.30  1  0    4    4

# unique函数一步搞定
unique(data3)
                  cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            6 167.6 123 3.92 3.440 18.30  1  0    4    4

数据翻转、更改、排序

# 数据框转置
t(mtcars[1:3,1:3])
     Mazda RX4 Mazda RX4 Wag Datsun 710
mpg         21            21       22.8
cyl          6             6        4.0
disp       160           160      108.0

# 翻转向量
letters
 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s" "t" "u" "v" "w" "x"
[25] "y" "z"
rev(letters)
 [1] "z" "y" "x" "w" "v" "u" "t" "s" "r" "q" "p" "o" "n" "m" "l" "k" "j" "i" "h" "g" "f" "e" "d" "c"
[25] "b" "a"

# 翻转数据框
women
   height weight
1      58    115
2      59    117
3      60    120
...
rev(row.names(women)) # 翻转行名后索引即可
 [1] "15" "14" "13" "12" "11" "10" "9"  "8"  "7"  "6"  "5"  "4"  "3"  "2"  "1" 
women[rev(row.names(women)), ]
   height weight
15     72    164
14     71    159
13     70    154
...

# 更改数据框某列的值
women$height*2.54
 [1] 147.32 149.86 152.40 154.94 157.48 160.02 162.56 165.10 167.64 170.18 172.72 175.26 177.80
[14] 180.34 182.88
data.frame(weight=women$height*2.54, women$weight)
   weight women.weight
1  147.32          115
2  149.86          117
3  152.40          120
...

# transform函数更高效，也可新增一列
transform(women, height=height*2.54)
   height weight
1  147.32    115
2  149.86    117
3  152.40    120
...

transform(women, cm=height*2.54)
   height weight     cm
1      58    115 147.32
2      59    117 149.86
3      60    120 152.40
...

# 数据排序
rivers
  [1]  735  320  325  392  524  450 1459  135  465  600  330  336  280  315  870  906  202  329  290
 [20] 1000  600  505 1450  840 1243  890  350  407  286  280  525  720  390  250  327  230  265  850
 [39]  210  630  260  230  360  730  600  306  390  420  291  710  340  217  281  352  259  250  470
 [58]  680  570  350  300  560  900  625  332 2348 1171 3710 2315 2533  780  280  410  460  260  255
 [77]  431  350  760  618  338  981 1306  500  696  605  250  411 1054  735  233  435  490  310  460
 [96]  383  375 1270  545  445 1885  380  300  380  377  425  276  210  800  420  350  360  538 1100
[115] 1205  314  237  610  360  540 1038  424  310  300  444  301  268  620  215  652  900  525  246
[134]  360  529  500  720  270  430  671 1770
rev(rivers)
  [1] 1770  671  430  270  720  500  529  360  246  525  900  652  215  620  268  301  444  300  310
 [20]  424 1038  540  360  610  237  314 1205 1100  538  360  350  420  800  210  276  425  377  380
 [39]  300  380 1885  445  545 1270  375  383  460  310  490  435  233  735 1054  411  250  605  696
 [58]  500 1306  981  338  618  760  350  431  255  260  460  410  280  780 2533 2315 3710 1171 2348
 [77]  332  625  900  560  300  350  570  680  470  250  259  352  281  217  340  710  291  420  390
 [96]  306  600  730  360  230  260  630  210  850  265  230  327  250  390  720  525  280  286  407
[115]  350  890 1243  840 1450  505  600 1000  290  329  202  906  870  315  280  336  330  600  465
[134]  135 1459  450  524  392  325  320  735

# 从小到大排序
sort(rivers)
  [1]  135  202  210  210  215  217  230  230  233  237  246  250  250  250  255  259  260  260  265
 [20]  268  270  276  280  280  280  281  286  290  291  300  300  300  301  306  310  310  314  315
 [39]  320  325  327  329  330  332  336  338  340  350  350  350  350  352  360  360  360  360  375
 [58]  377  380  380  383  390  390  392  407  410  411  420  420  424  425  430  431  435  444  445
 [77]  450  460  460  465  470  490  500  500  505  524  525  525  529  538  540  545  560  570  600
 [96]  600  600  605  610  618  620  625  630  652  671  680  696  710  720  720  730  735  735  760
[115]  780  800  840  850  870  890  900  900  906  981 1000 1038 1054 1100 1171 1205 1243 1270 1306
[134] 1450 1459 1770 1885 2315 2348 2533 3710
sort(rivers, decreasing = T)
  [1] 3710 2533 2348 2315 1885 1770 1459 1450 1306 1270 1243 1205 1171 1100 1054 1038 1000  981  906
 [20]  900  900  890  870  850  840  800  780  760  735  735  730  720  720  710  696  680  671  652
 [39]  630  625  620  618  610  605  600  600  600  570  560  545  540  538  529  525  525  524  505
 [58]  500  500  490  470  465  460  460  450  445  444  435  431  430  425  424  420  420  411  410
 [77]  407  392  390  390  383  380  380  377  375  360  360  360  360  352  350  350  350  350  340
 [96]  338  336  332  330  329  327  325  320  315  314  310  310  306  301  300  300  300  291  290
[115]  286  281  280  280  280  276  270  268  265  260  260  259  255  250  250  250  246  237  233
[134]  230  230  217  215  210  210  202  135

# 无法直接排序数据框，先对某一列排序后索引即可
sort(mtcars)
Error in `[.data.frame`(x, order(x, na.last = na.last, decreasing = decreasing)) : 
  选择了未定义的列
In addition: Warning message:
In xtfrm.data.frame(x) : cannot xtfrm data frames

# 根据行名排序
mtcars[sort(row.names(mtcars)), ]
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
...

rivers
  [1]  735  320  325  392  524  450 1459  135  465  600  330  336  280  315  870  906  202  329  290
 [20] 1000  600  505 1450  840 1243  890  350  407  286  280  525  720  390  250  327  230  265  850
 [39]  210  630  260  230  360  730  600  306  390  420  291  710  340  217  281  352  259  250  470
...
sort(rivers)
  [1]  135  202  210  210  215  217  230  230  233  237  246  250  250  250  255  259  260  260  265
 [20]  268  270  276  280  280  280  281  286  290  291  300  300  300  301  306  310  310  314  315
 [39]  320  325  327  329  330  332  336  338  340  350  350  350  350  352  360  360  360  360  375
...

# order，从小到大排序之后，数据在未排序状态下的序号，可用于索引
order(rivers)
  [1]   8  17  39 108 129  52  36  42  91 117 133  34  56  87  76  55  41  75  37 127 138 107  13  30
 [25]  72  53  29  19  49  61 103 124 126  46  94 123 116  14   2   3  35  18  11  65  12  81  51  27
 [49]  60  78 111  54  43 112 119 134  97 105 102 104  96  33  47   4  28  73  88  48 110 122 106 139
...
rivers[8]
[1] 135
rivers[17]
[1] 202

# 从大到小
order(-rivers)
  [1]  68  70  66  69 101 141   7  23  83  98  25 115  67 114  89 121  20  82  16  63 131  26  15  38
 [25]  24 109  71  79   1  90  44  32 137  50  85  58 140 130  40  64 128  80 118  86  10  21  45  59
 [49]  62  99 120 113 135  31 132   5  22  84 136  93  57   9  74  95   6 100 125  92  77 139 106 122
rivers[68]
[1] 3710

# 排序后索引即可对数据框排序
sort(mtcars$mpg)
 [1] 10.4 10.4 13.3 14.3 14.7 15.0 15.2 15.2 15.5 15.8 16.4 17.3 17.8 18.1 18.7 19.2 19.2 19.7 21.0
[20] 21.0 21.4 21.4 21.5 22.8 22.8 24.4 26.0 27.3 30.4 30.4 32.4 33.9
order(mtcars$mpg)
 [1] 15 16 24  7 17 31 14 23 22 29 12 13 11  6  5 10 25 30  1  2  4 32 21  3  9  8 27 26 19 28 18 20
mtcars[order(mtcars$mpg),]
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
...

#  可以多个排序条件
mtcars[order(mtcars$mpg, mtcars$disp),]
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
...
order(-rivers)
  [1]  68  70  66  69 101 141   7  23  83  98  25 115  67 114  89 121  20  82  16  63 131  26  15  38
 [25]  24 109  71  79   1  90  44  32 137  50  85  58 140 130  40  64 128  80 118  86  10  21  45  59
 [49]  62  99 120 113 135  31 132   5  22  84 136  93  57   9  74  95   6 100 125  92  77 139 106 122

数据框内数据计算

wp <- as.data.frame(WorldPhones)
head(wp)
     N.Amer Europe Asia S.Amer Oceania Africa Mid.Amer
1951  45939  21574 2876   1815    1646     89      555
1956  60423  29990 4708   2568    2366   1411      733
1957  64721  32510 5230   2695    2526   1546      773
...

# 计算每行列的和及平均值
rowSums(wp)
  1951   1956   1957   1958   1959   1960   1961 
 74494 102199 110001 118399 124801 133709 141700 
colMeans(wp)
    N.Amer     Europe       Asia     S.Amer    Oceania     Africa   Mid.Amer 
66747.5714 34343.4286  6229.2857  2772.2857  2625.0000  1484.0000   841.7143 
cbind(wp, Total=rowSums(wp))
     N.Amer Europe Asia S.Amer Oceania Africa Mid.Amer  Total
1951  45939  21574 2876   1815    1646     89      555  74494
1956  60423  29990 4708   2568    2366   1411      733 102199
1957  64721  32510 5230   2695    2526   1546      773 110001
...
new <- cbind(wp, Total=rowSums(wp))
rbind(new, Mean=colMeans(new))
       N.Amer   Europe     Asia   S.Amer Oceania Africa  Mid.Amer    Total
1951 45939.00 21574.00 2876.000 1815.000    1646     89  555.0000  74494.0
1956 60423.00 29990.00 4708.000 2568.000    2366   1411  733.0000 102199.0
1957 64721.00 32510.00 5230.000 2695.000    2526   1546  773.0000 110001.0
1958 68484.00 35218.00 6662.000 2845.000    2691   1663  836.0000 118399.0
1959 71799.00 37598.00 6856.000 3000.000    2868   1769  911.0000 124801.0
1960 76036.00 40341.00 8220.000 3145.000    3054   1905 1008.0000 133709.0
1961 79831.00 43173.00 9053.000 3338.000    3224   2005 1076.0000 141700.0
Mean 66747.57 34343.43 6229.286 2772.286    2625   1484  841.7143 115043.3

# 将一个FUN应用给数据框，包括自己编写的函数
apply(wp, 1, FUN = sum)
  1951   1956   1957   1958   1959   1960   1961 
 74494 102199 110001 118399 124801 133709 141700 
apply(wp, 2, sum)
  N.Amer   Europe     Asia   S.Amer  Oceania   Africa Mid.Amer 
  467233   240404    43605    19406    18375    10388     5892 
class(apply(wp, 2, sum))
[1] "numeric"
apply(wp, 1, sd)
    1951     1956     1957     1958     1959     1960     1961 
17309.22 22712.46 24362.16 25790.60 27116.58 28712.04 30213.56 

# 对列表操作，返回列表
lapply(state.center, length)
$x
[1] 50

$y
[1] 50
# 对列表操作，返回向量
sapply(state.center, length)
 x  y 
50 50 

# 第二个参数为因子
tapply(state.name, state.division, length)
       New England    Middle Atlantic     South Atlantic East South Central West South Central 
                 6                  3                  8                  4                  4 
East North Central West North Central           Mountain            Pacific 
                 5                  7                  8                  5 

# 数据中心化，各数据减去均值
x <- c(1,2,3,4,5)
x-mean(x)
[1] -2 -1  0  1  2

# 数据标准化，减去均值后再除以标准差
sd(x)
[1] 1.581139
(x-mean(x))/sd(x)
[1] -1.2649111 -0.6324555  0.0000000  0.6324555  1.2649111

scale(x, center = T, scale = F)
     [,1]
[1,]   -2
[2,]   -1
[3,]    0
[4,]    1
[5,]    2
attr(,"scaled:center")
[1] 3
scale(x, center = T, scale = T)
           [,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

heatmap(scale(state.x77, T,T))

`reshape2`包应用

`tidyr`包应用

作用：整理数据，合并数据，长宽数据互转

gather()：宽数据→长数据（建议使用更新的pivot_longer()和pivot_wider()）

spread()：长数据→宽数据

unite()：多列合并为一列

seperate()：一列拆分为多列

# gather()
td <- mtcars[1:2,1:3]
td
                   mpg cyl  disp
Mazda RX4         21.0   6 160.0
Mazda RX4 Wag     21.0   6 160.0
td <- data.frame(car=rownames(td),td)
td
                                car  mpg cyl  disp
Mazda RX4                 Mazda RX4 21.0   6 160.0
Mazda RX4 Wag         Mazda RX4 Wag 21.0   6 160.0

gather(td, key = 'Type', value = 'Value', mpg, cyl, disp)
            car Type Value
1     Mazda RX4  mpg    21
2 Mazda RX4 Wag  mpg    21
3     Mazda RX4  cyl     6
4 Mazda RX4 Wag  cyl     6
5     Mazda RX4 disp   160
6 Mazda RX4 Wag disp   160

# 选择转换的不同方法
gather(td, key = 'Type', value = 'Value', mpg:disp)
            car Type Value
1     Mazda RX4  mpg    21
2 Mazda RX4 Wag  mpg    21
3     Mazda RX4  cyl     6
4 Mazda RX4 Wag  cyl     6
5     Mazda RX4 disp   160
6 Mazda RX4 Wag disp   160

gather(td, key = 'Type', value = 'Value', mpg:disp, -cyl)
            car cyl Type Value
1     Mazda RX4   6  mpg    21
2 Mazda RX4 Wag   6  mpg    21
3     Mazda RX4   6 disp   160
4 Mazda RX4 Wag   6 disp   160

gather(td, key = 'Type', value = 'Value', 3:4)
            car mpg Type Value
1     Mazda RX4  21  cyl     6
2 Mazda RX4 Wag  21  cyl     6
3     Mazda RX4  21 disp   160
4 Mazda RX4 Wag  21 disp   160


# spread()
gtd <- gather(td, key = 'Type', value = 'Value', mpg:disp, -cyl)
gtd
            car cyl Type Value
1     Mazda RX4   6  mpg    21
2 Mazda RX4 Wag   6  mpg    21
3     Mazda RX4   6 disp   160
4 Mazda RX4 Wag   6 disp   160
spread(gtd, key = 'Type', value = 'Value')
            car cyl disp mpg
1     Mazda RX4   6  160  21
2 Mazda RX4 Wag   6  160  21

# separate()
x
     X
1 <NA>
2  1-2
3  2-3
4  3-4

separate(x, col = X, into = c('A','B'))
     A    B
1 <NA<NA>
2    1    2
3    2    3
4    3    4
separate(x, col = X, into = c('A','B'), sep = '-')
     A    B
1 <NA<NA>
2    1    2
3    2    3
4    3    4

# unite()
x
     A    B
1 <NA<NA>
2    1    2
3    2    3
4    3    4

unite(x, col=X)
      X
1 NA_NA
2   1_2
3   2_3
4   3_4
unite(x, col=X, sep = '-')
      X
1 NA-NA
2   1-2
3   2-3
4   3-4

# pivot_longer() 宽数据转换为长数据

`dplyr`包应用

# filter()过滤
> dplyr::filter(mtcars, mpg>30)
                mpg cyl disp  hp drat    wt  qsec vs am gear carb          model
Fiat 128       32.4   4 78.7  66 4.08 2.200 19.47  1  1    4    1       Fiat 128
Honda Civic    30.4   4 75.7  52 4.93 1.615 18.52  1  1    4    2    Honda Civic
Toyota Corolla 33.9   4 71.1  65 4.22 1.835 19.90  1  1    4    1 Toyota Corolla
Lotus Europa   30.4   4 95.1 113 3.77 1.513 16.90  1  1    5    2   Lotus Europa
> dplyr::filter(mtcars, mpg>30 & wt>2)
          mpg cyl disp hp drat  wt  qsec vs am gear carb    model
Fiat 128 32.4   4 78.7 66 4.08 2.2 19.47  1  1    4    1 Fiat 128

# distinct()去除重复
> dplyr::distinct(rbind(mtcars[1:2,], mtcars[1:3,]))
               mpg cyl disp  hp drat    wt  qsec vs am gear carb         model
Mazda RX4     21.0   6  160 110 3.90 2.620 16.46  0  1    4    4     Mazda RX4
Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4 Mazda RX4 Wag
Datsun 710    22.8   4  108  93 3.85 2.320 18.61  1  1    4    1    Datsun 710

# slice()切片
> dplyr::slice(mtcars, 1:3)
               mpg cyl disp  hp drat    wt  qsec vs am gear carb         model
Mazda RX4     21.0   6  160 110 3.90 2.620 16.46  0  1    4    4     Mazda RX4
Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4 Mazda RX4 Wag
Datsun 710    22.8   4  108  93 3.85 2.320 18.61  1  1    4    1    Datsun 710

# sample_n()抽样，sample_frac()按比例抽样
> dplyr::sample_n(mtcars, 3)
               mpg cyl  disp  hp drat    wt  qsec vs am gear carb         model
Porsche 914-2 26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2 Porsche 914-2
Fiat 128      32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1      Fiat 128
Mazda RX4 Wag 21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4 Mazda RX4 Wag

> length(row.names(mtcars))
[1] 32
> dplyr::sample_frac(mtcars, 0.1)
              mpg cyl  disp  hp drat    wt qsec vs am gear carb        model
Ferrari Dino 19.7   6 145.0 175 3.62 2.770 15.5  0  1    5    6 Ferrari Dino
AMC Javelin  15.2   8 304.0 150 3.15 3.435 17.3  0  0    3    2  AMC Javelin
Merc 450SE   16.4   8 275.8 180 3.07 4.070 17.4  0  0    3    3   Merc 450SE

# arrange()排序
> dplyr::arrange(mtcars, mpg)
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb               model
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4  Cadillac Fleetwood
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4 Lincoln Continental
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4          Camaro Z28
...
> dplyr::arrange(mtcars, -mpg)
> dplyr::arrange(mtcars, desc(mpg))
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb               model
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1      Toyota Corolla
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1            Fiat 128
Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2         Honda Civic
...

# select()数据框取子集

# summarise()描述性统计，还有summarise_all(), summarise_at(), summarise_each()等
> summarise(mtcars, avg=mean(mpg))
       avg
1 20.09062
> summarise(mtcars, avg=mean(cyl))
     avg
1 6.1875
> summarise(mtcars, sum=sum(mpg))
    sum
1 642.9

# % > %链式操作符，管道符，前一个运算结果作为后一个的输入
> mtcars %>% summarise(sum=sum(mpg))
    sum
1 642.9
> head(mtcars, 10) %>% tail(3)
           mpg cyl  disp  hp drat   wt qsec vs am gear carb
Merc 240D 24.4   4 146.7  62 3.69 3.19 20.0  1  0    4    2
Merc 230  22.8   4 140.8  95 3.92 3.15 22.9  1  0    4    2
Merc 280  19.2   6 167.6 123 3.92 3.44 18.3  1  0    4    4

# group_by()分组
> dplyr::group_by(mtcars, cyl)
# A tibble: 32 x 11
# Groups:   cyl [3]
     mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1  21       6  160    110  3.9   2.62  16.5     0     1     4     4
 2  21       6  160    110  3.9   2.88  17.0     0     1     4     4
 3  22.8     4  108     93  3.85  2.32  18.6     1     1     4     1
 4  21.4     6  258    110  3.08  3.22  19.4     1     0     3     1
 5  18.7     8  360    175  3.15  3.44  17.0     0     0     3     2
 6  18.1     6  225    105  2.76  3.46  20.2     1     0     3     1
 7  14.3     8  360    245  3.21  3.57  15.8     0     0     3     4
 8  24.4     4  147.    62  3.69  3.19  20       1     0     4     2
 9  22.8     4  141.    95  3.92  3.15  22.9     1     0     4     2
10  19.2     6  168.   123  3.92  3.44  18.3     1     0     4     4
# ... with 22 more rows

> dplyr::group_by(mtcars, cyl) %>% summarise()
# A tibble: 3 x 1
    cyl
  <dbl>
1     4
2     6
3     8
> dplyr::group_by(mtcars, cyl) %>% summarise(ave=mean(mpg))
# A tibble: 3 x 2
    cyl   ave
  <dbl> <dbl>
1     4  26.7
2     6  19.7
3     8  15.1
> dplyr::group_by(mtcars, cyl) %>% summarise(ave=mean(mpg)) %>% arrange(ave)
# A tibble: 3 x 2
    cyl   ave
  <dbl> <dbl>
1     8  15.1
2     6  19.7
3     4  26.7

# mutate()加入新变量
> mt <- mutate(mtcars, new = mpg+disp)
> head(mt,3)
               mpg cyl disp  hp drat    wt  qsec vs am gear carb   new
Mazda RX4     21.0   6  160 110 3.90 2.620 16.46  0  1    4    4 181.0
Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4 181.0
Datsun 710    22.8   4  108  93 3.85 2.320 18.61  1  1    4    1 130.8

# 整合x, y两个表格，类似于集合的操作
> a
  x1 x2
1  A  1
2  B  2
3  C  3
> b
  x1    x3
1  A  TRUE
2  B FALSE
3  D  TRUE

# 左右链接，以某列为标准
> dplyr::left_join(a,b,by="x1")
  x1 x2    x3
1  A  1  TRUE
2  B  2 FALSE
3  C  3    NA
> dplyr::right_join(a, b, 'x1')
  x1 x2    x3
1  A  1  TRUE
2  B  2 FALSE
3  D NA  TRUE

# 内链接取x1交集，外链接取全集
> inner_join(a, b, 'x1')
  x1 x2    x3
1  A  1  TRUE
2  B  2 FALSE
> dplyr::full_join(a,b,by="x1")
  x1 x2    x3
1  A  1  TRUE
2  B  2 FALSE
3  C  3    NA
4  D NA  TRUE

# 半链接取交集，反链接取补集
> dplyr::semi_join(a,b,by="x1")
  x1 x2
1  A  1
2  B  2
> dplyr::anti_join(a,b,by="x1")
  x1 x2
1  C  3

# 数据框交集、补集、并集操作
> car1
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
> car2
                   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2

# 取交集
> intersect(car1, car2)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2

# 取并集，不包括重复项
> union(car1, car2)
                   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2

# 取并集，包括重复项
> union_all(car1, car2)
                       mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4             21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag         21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710...3        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive...4    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout...5 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Datsun 710...6        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive...7    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout...8 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant               18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360            14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D             24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2

# 取car1补集
> setdiff(car1, car2)
              mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4      21   6  160 110  3.9 2.620 16.46  0  1    4    4
Mazda RX4 Wag  21   6  160 110  3.9 2.875 17.02  0  1    4    4
# 取car2补集
> setdiff(car2, car1)
            mpg cyl  disp  hp drat   wt  qsec vs am gear carb
Valiant    18.1   6 225.0 105 2.76 3.46 20.22  1  0    3    1
Duster 360 14.3   8 360.0 245 3.21 3.57 15.84  0  0    3    4
Merc 240D  24.4   4 146.7  62 3.69 3.19 20.00  1  0    4    2

R函数

函数的选项参数

函数选项类型

输入控制
输出控制
调节

数学统计函数

d概率密度函数

p分布函数

q分布函数的反函数

r随机数函数

如正态分布函数norm()加上这4个前缀：

dnorm(), pnorm(), qnorm(), rnorm()

# 生成10个符合正态分布的随机数
> rnorm(n=10, mean=5, sd=1)
 [1] 6.329543 4.232249 6.033233 4.837938 3.984797 5.299808 5.202550 2.070615 5.830317 5.936402
> round(rnorm(n=10, mean=5, sd=1))
 [1] 5 4 5 5 5 4 4 5 6 3

# runif()生成随机数
> runif(5)
[1] 0.04366431 0.34873894 0.29974298 0.39224632 0.08980883
> runif(5, 0, 10)
[1] 0.83081773 2.22285093 0.07499167 4.59838033 6.37444314
> runif(5, min=0, max=10)
[1] 4.5102458 0.3494494 9.4862655 4.6722139 0.5983668
> round(runif(5, min=0, max=10))
[1]  1  5  8  6 10

# 固定生成的随机数
> set.seed(10)
> round(runif(5, min=0, max=10))
[1] 5 3 4 7 1
> round(runif(5, min=0, max=10))
[1] 2 3 3 6 4
> round(runif(5, min=0, max=10))
[1] 7 6 1 6 4
> set.seed(10)
> round(runif(5, min=0, max=10))
[1] 5 3 4 7 1

描述性统计函数

频数统计函数

# table()直接统计频数
> table(mtcars$cyl)

 4  6  8 
11  7 14 
# prop.table()计算频率
> prop.table(table(mtcars$cyl))

      4       6       8 
0.34375 0.21875 0.43750 

# split()将数据分组(list)
> split(mtcars, mtcars$cyl)
$`4`
                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Datsun 710     22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2

$`6`
                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4      21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4

$`8`
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
> split(mtcars, mtcars$cyl)['4']
$`4`
                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Datsun 710     22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2

# 计算连续值频数和频率
> range(mtcars$mpg)
[1] 10.4 33.9
> seq(10, 40, 5)
[1] 10 15 20 25 30 35 40
> cut(mtcars$mpg, seq(10, 35, 5))
 [1] (20,25] (20,25] (20,25] (20,25] (15,20] (15,20] (10,15] (20,25] (20,25] (15,20] (15,20] (15,20]
[13] (15,20] (15,20] (10,15] (10,15] (10,15] (30,35] (30,35] (30,35] (20,25] (15,20] (15,20] (10,15]
[25] (15,20] (25,30] (25,30] (30,35] (15,20] (15,20] (10,15] (20,25]
Levels: (10,15] (15,20] (20,25] (25,30] (30,35]

> table(cut(mtcars$mpg, seq(10, 35, 5)))

(10,15] (15,20] (20,25] (25,30] (30,35] 
      6      12       8       2       4 

> prop.table(table(cut(mtcars$mpg, seq(10, 35, 5))))

(10,15] (15,20] (20,25] (25,30] (30,35] 
 0.1875  0.3750  0.2500  0.0625  0.1250 

# 二维，talbe() xtabs()
> library(vcd) 
> Arthritis
   ID Treatment    Sex Age Improved
1  57   Treated   Male  27     Some
2  46   Treated   Male  29     None
3  77   Treated   Male  30     None
                                 
> table(Arthritis$Treatment, Arthritis$Improved)
         
          None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
                                 
> xtabs(~ Treatment + Improved, data = Arthritis)
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

> with(data = Arthritis, {table(Treatment, Improved)})
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
 
# 边际频数及统计，1表示按行统计，2表示按列统计
# margin.table() addmargins()
> x <- xtabs(~ Treatment + Improved, data = Arthritis)
> margin.table(x)
[1] 84
> margin.table(x, 1)
Treatment
Placebo Treated 
     43      41 
> margin.table(x, 2)
Improved
  None   Some Marked 
    42     14     28 
                                 
> prop.table(x, 1)
         Improved
Treatment      None      Some    Marked
  Placebo 0.6744186 0.1627907 0.1627907
  Treated 0.3170732 0.1707317 0.5121951
> prop.table(x, 2)
         Improved
Treatment      None      Some    Marked
  Placebo 0.6904762 0.5000000 0.2500000
  Treated 0.3095238 0.5000000 0.7500000
> prop.table(x)
         Improved
Treatment       None       Some     Marked
  Placebo 0.34523810 0.08333333 0.08333333
  Treated 0.15476190 0.08333333 0.25000000

> addmargins(x)
         Improved
Treatment None Some Marked Sum
  Placebo   29    7      7  43
  Treated   13    7     21  41
  Sum       42   14     28  84
> addmargins(x, 1)
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
  Sum       42   14     28
> addmargins(x, 2)
         Improved
Treatment None Some Marked Sum
  Placebo   29    7      7  43
  Treated   13    7     21  41

# 三维ftable()
> y <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
> y
, , Sex = Female

         Improved
Treatment None Some Marked
  Placebo   19    7      6
  Treated    6    5     16

, , Sex = Male

         Improved
Treatment None Some Marked
  Placebo   10    0      1
  Treated    7    2      5

> ftable(y)
                   Sex Female Male
Treatment Improved                
Placebo   None             19   10
          Some              7    0
          Marked            6    1
Treated   None              6    7
          Some              5    2
          Marked           16    5

独立性检验函数

独立性检验是根据频数信息判断两类因子彼此相关或者独立的假设检验。

独立性检验算法：

卡方检验
Fisher检验
Cochran-Mantel-Haenszel检验

假设检验（Hypothesis Testing）：由抽样判断总体。

原假设：没有发生；

备择假设：发生了

p value：原假设为真时，得到最大的或者超出所得到的的检验统计计算量的概率

p值一般定到0.05，当p<0.05时，拒绝原假设，p>0.05时，不拒绝原假设（根据情况调整p值）

> Arthritis
   ID Treatment    Sex Age Improved
1  57   Treated   Male  27     Some
2  46   Treated   Male  29     None
3  77   Treated   Male  30     None

> table(Arthritis$Treatment, Arthritis$Improved)
         
          None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
> mytable <- table(Arthritis$Treatment, Arthritis$Improved)

# 卡方检验，p<0.05，拒绝原假设，Treatment和Improved有关
> chisq.test(mytable)

	Pearson's Chi-squared test

data:  mytable
X-squared = 13.055, df = 2, p-value = 0.001463

# p>0.05，不拒绝原假设，Sex和Improved独立
> mytable <- table(Arthritis$Sex, Arthritis$Improved)
> chisq.test(mytable)

	Pearson's Chi-squared test

data:  mytable
X-squared = 4.8407, df = 2, p-value = 0.08889

Warning message:
In chisq.test(mytable) : Chi-squared approximation may be incorrect

# Fisher检验，p<0.05，拒绝原假设，Treatment和Improved有关
> mytable <- table(Arthritis$Treatment, Arthritis$Improved)
> fisher.test(mytable)

	Fisher's Exact Test for Count Data

data:  mytable
p-value = 0.001393
alternative hypothesis: two.sided

# Cochran-Mantel-Haenszel检验，p<0.05，拒绝原假设，Treatment和Improved有关
> mytable <- xtabs(~ Treatment+Improved+Sex, data = Arthritis)
> mytable
, , Sex = Female

         Improved
Treatment None Some Marked
  Placebo   19    7      6
  Treated    6    5     16

, , Sex = Male

         Improved
Treatment None Some Marked
  Placebo   10    0      1
  Treated    7    2      5

> mantelhaen.test(mytable)

	Cochran-Mantel-Haenszel test

data:  mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647

# 顺序很重要
> mytable <- xtabs(~ Treatment+Sex+Improved, data = Arthritis)
> mantelhaen.test(mytable)

	Mantel-Haenszel chi-squared test with continuity correction

data:  mytable
Mantel-Haenszel X-squared = 2.0863, df = 1, p-value = 0.1486
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.8566711 8.0070521
sample estimates:
common odds ratio 
         2.619048 

> mytable <- xtabs(~ Sex+Improved+Treatment, data = Arthritis)
> mantelhaen.test(mytable)

	Cochran-Mantel-Haenszel test

data:  mytable
Cochran-Mantel-Haenszel M^2 = 6.846, df = 2, p-value = 0.03261

绘图函数

见第三章

自定义函数

# 自定义计算偏度和丰度的函数
myfun <- function(x, na.omit=FALSE){
  if(na.omit)
    x <- x[!is.na(x)]
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  skew <- sum(x-m^3/s^3)/n
  kurt <- sum(x-m^4/s^4)/n-3
  return(c(n=m, mean=m, stdev=s, skew=skew, kurtosis=kurt))
}

# 循环
> for(i in 1:10){print('11')}

> i=1
> while(i<=10){print('2');i=i+1}

> ifelse(i>10, print('yes'), print('No'))

数据分析实战

线性回归

回归（Regression）：用用预测变量（自变量、解释变量）来预测响应变量（因变量、结果变量）

回归	用途
简单线性	用一个量化的解释变量预测一个量化的响应变量
多项式	用一个量化的解释变量预测一个量化的响应变量，模型的关系是n阶多项式
多元线性	用两个或多个量化的解释变量预测一个量化的响应变量
多变量	用一个或多个解释变量预测多个响应变量
Logistic	用一个或多个解释变量预测一个类别型响应变量
泊松	用一个或多个解释变量预测一个代表频数的响应变量
Cox比例风险	用一个或多个解释变量预测一个事件 ( 死亡、失败或旧病复发)发生的时间
时间序列	对误差项相关的时间序列数据建模
非线性	用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
非参数	用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式,不事先设定
稳健	用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

线性回归：

简单线性回归
多项式
多元线性
多变量
Logistic

简单线性回归

普通最小二乘回归法（OLS）

lm()函数进行线性回归

> women
   height weight
1      58    115
2      59    117
3      60    120

> fit <- lm(weight ~ height, data = women)
> fit

Call:
lm(formula = weight ~ height, data = women)

Coefficients:
(Intercept)       height  
     -87.52         3.45 

# 回归分析结果
> summary(fit)

Call:
lm(formula = weight ~ height, data = women) # 回归所用公式

Residuals: # 残差：真实值与预测值的差，越小越精确，残差分布
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)     # p value
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 *** # 截距
height        3.45000    0.09114   37.85 1.09e-14 *** # 系数
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.525 on 13 degrees of freedom # 残差的标准误差，越小越好
Multiple R-squared:  0.991,	Adjusted R-squared:  0.9903 # R2判定系数，回归模型能解释的响应变量的方差比例
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14  # F统计量，说明模型是否显著
# 一般先看F统计量，p>0.05模型就没有价值，需重新拟合
# 再看R2，看模型能够解释多少变量

响应变量~预测变量1+预测变量2+预测变量3

拟合线性模型的常用函数

> fit

Call:
lm(formula = weight ~ height, data = women)

Coefficients:
(Intercept)       height  
     -87.52         3.45  

# 列出拟合模型的截距和斜率
> coefficients(fit)
(Intercept)      height 
  -87.51667     3.45000 

# 给出置信区间
> confint(fit)
                  2.5 %     97.5 %
(Intercept) -100.342655 -74.690679
height         3.253112   3.646888
> confint(fit, level = 0.5)
                  25 %       75 %
(Intercept) -91.635892 -83.397441
height        3.386767   3.513233

# 给出预测值
> fitted(fit)
       1        2        3        4        5        6        7        8        9       10 
112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 140.1833 143.6333 
      11       12       13       14       15 
147.0833 150.5333 153.9833 157.4333 160.8833 

# 给出残差
> residuals(fit)
          1           2           3           4           5           6           7           8 
 2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 -1.28333333 -1.73333333 
          9          10          11          12          13          14          15 
-1.18333333 -1.63333333 -1.08333333 -0.53333333  0.01666667  1.56666667  3.11666667 
 
> women$weight-fitted(fit)
          1           2           3           4           5           6           7           8 
 2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 -1.28333333 -1.73333333 
          9          10          11          12          13          14          15 
-1.18333333 -1.63333333 -1.08333333 -0.53333333  0.01666667  1.56666667  3.11666667 

> women1 <- women
# 利用拟合模型对新数据集的响应变量进行预测
> predict(fit, women1)
       1        2        3        4        5        6        7        8        9       10 
112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 140.1833 143.6333 
      11       12       13       14       15 
147.0833 150.5333 153.9833 157.4333 160.8833 

# 生成评价拟合模型的诊断图
> plot(fit)

多项式回归

# 多项式用I()括起来
> fit2 <- lm(weight~height+I(height^2), data=women)
> fit2

Call:
lm(formula = weight ~ height + I(height^2), data = women)

Coefficients: # 给出各项系数
(Intercept)       height  I(height^2)  
  261.87818     -7.34832      0.08306  

# 拟合效果优于简单线性回归
> summary(fit2)

Call:
lm(formula = weight ~ height + I(height^2), data = women)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50941 -0.29611 -0.00941  0.28615  0.59706 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
height       -7.34832    0.77769  -9.449 6.58e-07 ***
I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3841 on 12 degrees of freedom
Multiple R-squared:  0.9995,	Adjusted R-squared:  0.9994 
F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16

> fit3 <- lm(weight~height+I(height^2)+I(height^3), data=women)
> summary(fit3)

Call:
lm(formula = weight ~ height + I(height^2) + I(height^3), data = women)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.40677 -0.17391  0.03091  0.12051  0.42191 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -8.967e+02  2.946e+02  -3.044  0.01116 * 
height       4.641e+01  1.366e+01   3.399  0.00594 **
I(height^2) -7.462e-01  2.105e-01  -3.544  0.00460 **
I(height^3)  4.253e-03  1.079e-03   3.940  0.00231 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2583 on 11 degrees of freedom
Multiple R-squared:  0.9998,	Adjusted R-squared:  0.9997 
F-statistic: 1.679e+04 on 3 and 11 DF,  p-value: < 2.2e-16

# 查看拟合效果
> plot(women)
> abline(fit)
> lines(women$height, fitted(fit2), col='red')
> lines(women$height, fitted(fit3), col='blue')

多元线性回归

回归系数的含义：其他预测变量保持不变时，一个预测变量增加一个单位后，因变量的增量。

> state <- as.data.frame(state.x77[,c('Murder', 'Life Exp', 'Illiteracy', 'Income', 'HS Grad')])
> colnames(state) <- c('Murder', 'Life_Exp', 'Illiteracy', 'Income', 'HS_Grad')
 
# 加入多个变量
> fit <- lm(Murder ~ Life_Exp+Illiteracy+Income+HS_Grad, data=state)
> fit

Call:
lm(formula = Murder ~ Life_Exp + Illiteracy + Income + HS_Grad, 
    data = state)

Coefficients:
(Intercept)     Life_Exp   Illiteracy       Income      HS_Grad  
 116.507133    -1.662190     2.786755     0.000718     0.042163  

> summary(fit)

Call:
lm(formula = Murder ~ Life_Exp + Illiteracy + Income + HS_Grad, 
    data = state)

Residuals:
   Min     1Q Median     3Q    Max 
-3.715 -1.524  0.113  1.737  3.848 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.165e+02  1.976e+01   5.897 4.44e-07 ***
Life_Exp    -1.662e+00  2.822e-01  -5.891 4.53e-07 ***
Illiteracy   2.787e+00  6.708e-01   4.154 0.000144 ***
Income       7.180e-04  6.024e-04   1.192 0.239535    
HS_Grad      4.216e-02  5.736e-02   0.735 0.466095    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.028 on 45 degrees of freedom
Multiple R-squared:  0.7229,	Adjusted R-squared:  0.6983 
F-statistic: 29.35 on 4 and 45 DF,  p-value: 4.961e-12

# 查看系数
> coef(fit)
  (Intercept)      Life_Exp    Illiteracy        Income       HS_Grad 
 1.165071e+02 -1.662190e+00  2.786755e+00  7.179837e-04  4.216343e-02 

# 回归模型比较
> fit1 <- lm(Murder ~ Life_Exp+Illiteracy+Income+HS_Grad, data=state)
> fit2 <- lm(Murder ~ Life_Exp+Illiteracy, data=state)
> summary(fit2)

Call:
lm(formula = Murder ~ Life_Exp + Illiteracy, data = state)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3725 -1.5867  0.0284  1.1919  4.3641 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 114.2162    19.6698   5.807 5.27e-07 ***
Life_Exp     -1.5446     0.2716  -5.688 7.96e-07 ***
Illiteracy    2.2557     0.5981   3.772 0.000453 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.063 on 47 degrees of freedom
Multiple R-squared:  0.7004,	Adjusted R-squared:  0.6876 
F-statistic: 54.94 on 2 and 47 DF,  p-value: 4.998e-13

# AIC()函数比较拟合模型，越小越好，表明模型用较少的参数就可以得到很好的拟合度
> AIC(fit1, fit2)
     df      AIC
fit1  6 219.3144
fit2  4 219.2232

有交互项的多元线性回归

有时候变量之间不是独立的，有交互关系，如车辆马力与车重有关系.

若两个预测变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另一个预测变量的水平。

> fit <- lm(mpg~hp+wt+hp:wt, data=mtcars)
> summary(fit)

Call:
lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars) # 不清楚是什么交互关系，用:连接

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0632 -1.6491 -0.7362  1.4211  4.5513 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
hp          -0.12010    0.02470  -4.863 4.04e-05 ***
wt          -8.21662    1.26971  -6.471 5.20e-07 ***
hp:wt        0.02785    0.00742   3.753 0.000811 *** # hp与wt交互显著
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.153 on 28 degrees of freedom
Multiple R-squared:  0.8848,	Adjusted R-squared:  0.8724 
F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13

AIC()只能比较两个模型

如果变量很多，可以采用逐步回归法和全子集回归法

逐步回归法：模型会一次添加或删除变量，向前或向后逐步回归

全子集回归法：计算所有模型后选出最优模型，更全面

逐步回归法

> library(MASS)
> states <- as.data.frame(state.x77[,c("Murder", "Population",
+                                      "Illiteracy", "Income", "Frost")])
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,
+           data=states)
> stepAIC(fit, direction = "backward")
Start:  AIC=97.75
Murder ~ Population + Illiteracy + Income + Frost

             Df Sum of Sq    RSS     AIC
- Frost       1     0.021 289.19  95.753
- Income      1     0.057 289.22  95.759
<none>                    289.17  97.749
- Population  1    39.238 328.41 102.111
- Illiteracy  1   144.264 433.43 115.986

Step:  AIC=95.75
Murder ~ Population + Illiteracy + Income

             Df Sum of Sq    RSS     AIC
- Income      1     0.057 289.25  93.763
<none>                    289.19  95.753
- Population  1    43.658 332.85 100.783
- Illiteracy  1   236.196 525.38 123.605

Step:  AIC=93.76
Murder ~ Population + Illiteracy

             Df Sum of Sq    RSS     AIC
<none>                    289.25  93.763
- Population  1    48.517 337.76  99.516
- Illiteracy  1   299.646 588.89 127.311

Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)

Coefficients:
(Intercept)   Population   Illiteracy  
  1.6515497    0.0002242    4.0807366

全子集回归法

> library(leaps)
Warning message:

> states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
> leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
> plot(leaps, scale="adjr2")

回归诊断

模型是否最佳？是否经得起更多数据的检验？

满足OLS模型的统计假设才能够用lm()进行拟合

正态性：对于固定的自变量值，因变量成正态分布
独立性：因变量之间相互独立
线性：因变量和自变量之间为线性关系
同方差性：因变量的方差不随自变量的水平不同而变化。也称不变方差

plot(fit)的四幅图即是对这四个条件的反应

opar <- par(no.readonly=TRUE)
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)
par(opar)

残差——拟合图

ANOVA

Analysis of Variance，ANOVA，也称变异数分析，一种统计检验，用于分析两个以上样本之间的均值的差异。

方差分析分类：

单因素方差分析ANOVA（组内、组间）
双因素方差分析ANOVA
协方差分析ANCOVA
多元方差分析MANOVA
多元协方差分析MANCOVA

aov()和lm()均可进行方差分析，aov()用法类似

常见的公式

公式中各项的顺序很重要。

one-way ANOVA

探究单个单个因素（自变量）对响应变量的影响。

one-way ANOVA，单因素方差分析，一个自变量（independent variable）

适用情况：自变量至少有3个水平（level，即类型为因子factor），即至少3个不同组或类别，因变量（dependent variable）为数值型（numeric）

ANOVA可检测因变量是否随自变量的水平变化而变化，比如

自变量为微信使用频率，分为低、中、高三个组，探究其晚上睡眠时长是否有差异。
自变量为汽水品牌，收集可口可乐、百事、芬达、康师傅的数据，分析其每100 mL汽水的价格是否有差异。
自变量为肥料种类，使用3种肥料对作物施肥，分析作物产量是否有差异。

ANOVA的零假设（null hypothesis，H₀）：不同组之间的均值没有差异

备择假设（alternative hypothesis， H_a）：至少有一个组与因变量的总体均值显著不同

原理：计算各组的均值是否与因变量总体均值不同，如果有任何一个组的均值与总体均值显著不同，则拒绝零假设。ANOVA使用F检验（F-test）进行统计的显著性分析。由于误差是针对整个组比较计算的，所以F检验允许一次比较多个均值，而不是针对每个单独的双向比较（t-test）。F检验将每组均值的方差与总体组的方差进行比较，如果组内方差小于组间方差，则F检验得到的F值很大，因此观察到的差异是真实存在的可能性更高，而不是偶然现象。

ANOVA假设条件

观测之间相互独立：使用有效的统计学方法收集数据，观测之间没有隐藏的联系
响应变量符合正态分布：因变量的值符合正态分布
同方差性：被比较的每组的方差是相似的，如果组间方差不同，则方差分析可能不适合

实例

# 作物数据
# density：种植密度，1：低，2：高
# block：种植区域，1,2,3,4
# fertilizer：肥料种类：1,2,3
# yield：产量（bushels/acre，蒲式耳/英亩）

# 读入数据时，直接转换为 factor
> oneway <- read.csv('crop.data.csv', colClasses = c('factor', 'factor', 'factor', 'numeric'))
> head(oneway)
  density block fertilizer    yield
1       1     1          1 177.2287
2       2     2          1 177.5500
3       1     3          1 176.4085
4       2     4          1 177.7036
5       1     1          1 177.1255
6       2     2          1 176.7783
> summary(oneway)
 density block  fertilizer     yield      
 1:48    1:24   1:32       Min.   :175.4  
 2:48    2:24   2:32       1st Qu.:176.5  
         3:24   3:32       Median :177.1  
         4:24              Mean   :177.0  
                           3rd Qu.:177.4  
                           Max.   :179.1  

# 探究肥料与产量的关系
> fit <- aov(yield~fertilizer, data = oneway)
> summary(fit)
            Df Sum Sq Mean Sq F value Pr(>F)    
fertilizer   2   6.07  3.0340   7.863  7e-04 ***
Residuals   93  35.89  0.3859                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANOVA 输出内容解释

第一列：自变量和模型残差（model residuals），也即模型误差（model error）
Df：为自变量的自由度（自变量level-1）和残差的自由度（观测总数-1-？？？？），fertilizer自由度为3-1=2，残差为96-1-2=93
Sum Sq：由该变量解释的组均值和总体均值之间的平方和（sum of squares），fertilizer平方和为6.07，残差为35.89
Mean Sq：平方和的均值，即平方和÷自由度
F value：来自F检验（F-test）的检验统计量，每个自变量的平均方差÷残差的平均方差，F值越大，与自变量相关的变化越有可能是真实的，而不是偶然的
Pr(>F)：F统计量的p值，表示，如果组均值之间没有差异的零假设为真时，F值发生的可能性

由于自变量fertilizer的p < 0.05，显著，拒绝零假设，所以肥料种类对作物产量有显著影响。

ANOVA 能够分析自变量的水平之间是否存在差异，但哪些差异是显著的却不知道。

可以使用TukeyHSD()函数进行 Tukey’s Honestly-Significant Difference post-hoc test

Tukey test 对组进行两两比较，并使用保守误差（conservative error estimate）估计找出在统计上不同的组

> TukeyHSD(fit)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = yield ~ fertilizer, data = oneway)

$fertilizer
         diff         lwr       upr     p adj
2-1 0.1761687 -0.19371896 0.5460564 0.4954705
3-1 0.5991256  0.22923789 0.9690133 0.0006125
3-2 0.4229569  0.05306916 0.7928445 0.0208735

> aggregate(oneway$yield, by=list(oneway$fertilizer), FUN=mean)
  Group.1        x
1       1 176.7570
2       2 176.9332
3       3 177.3562

输出内容解释

Fit：检验的模型

diff：两组的均值差

lwr, upr：95%置信区间的上下限

p adj：调整后的p值

结果表明，3-2和3-1之间的差异显著（p < 0.05），肥料3的作物平均产量显著高于肥料1和2；

而肥料1和2之间的作物平均产量差异不显著。

ANOVA 总结方式

We found a statistically-significant difference in average crop yield according to fertilizer type (f(2)=9.073, p < 0.001). A Tukey post-hoc test revealed significant pairwise differences between fertilizer types 3 and 2, with an average difference of 0.42 bushels/acre (p < 0.05) and between fertilizer types 3 and 1, with an average difference of 0.59 bushels/acre (p < 0.01).

multcomp包中的cholesterol数据集

自变量：5种疗法trt

响应变量：问卷得分response

> library(multcomp)
> head(cholesterol)
    trt response
1 1time   3.8612
2 1time  10.3868
3 1time   5.9059
4 1time   3.0609
5 1time   7.7204
6 1time   2.7139

# 频数统计
> table(cholesterol$trt)

 1time 2times 4times  drugD  drugE 
    10     10     10     10     10 

# 分组统计均值和标准差
> aggregate(cholesterol$response, by=list(cholesterol$trt), FUN=mean)
  Group.1        x
1   1time  5.78197
2  2times  9.22497
3  4times 12.37478
4   drugD 15.36117
5   drugE 20.94752

> aggregate(cholesterol$response, by=list(cholesterol$trt), FUN=sd)
  Group.1        x
1   1time 2.878113
2  2times 3.483054
3  4times 2.923119
4   drugD 3.454636
5   drugE 3.345003

# 方差分析
> fit <- aov(response~trt, data = cholesterol)
> fit
Call:
   aov(formula = response ~ trt, data = cholesterol)

Terms:
                      trt Residuals
Sum of Squares  1351.3690  468.7504
Deg. of Freedom         4        45

Residual standard error: 3.227488
Estimated effects may be unbalanced

> summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 *** # F值越大，组间差异越大，p值越小，F值越可靠
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 检测数据是否满足假设条件
> plot(fit)


# lm()函数方差分析
> fit_lm <- lm(response~trt, data = cholesterol)
> summary(fit_lm)

Call:
lm(formula = response ~ trt, data = cholesterol)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5418 -1.9672 -0.0016  1.8901  6.6008 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.782      1.021   5.665 9.78e-07 ***
trt2times      3.443      1.443   2.385   0.0213 *  
trt4times      6.593      1.443   4.568 3.82e-05 ***
trtdrugD       9.579      1.443   6.637 3.53e-08 ***
trtdrugE      15.166      1.443  10.507 1.08e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.227 on 45 degrees of freedom
Multiple R-squared:  0.7425,	Adjusted R-squared:  0.7196 
F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13 # F值与p值和aov()函数得出的相同

单因素协方差分析

litter怀孕小鼠数据集

因变量：给药剂量dose，0, 5, 50, 500

响应变量：每窝幼崽平均出生体重weight

协变量：怀孕时间gesttime、产仔数

> library(multcomp)
> head(litter)
  dose weight gesttime number
1    0  28.05     22.5     15
2    0  33.33     22.5     14
3    0  36.37     22.0     14
4    0  35.52     22.0     13
5    0  36.77     21.5     15
6    0  29.60     23.0      5

# 频数统计
> table(litter$dose)

  0   5  50 500 
 20  19  18  17 
> table(litter$gesttime)

21.5   22 22.5   23 
  20   24   27    3 

# 分组统计平均值，差异不明显
> aggregate(litter$weight, by=list(litter$dose), FUN=mean)
  Group.1        x
1       0 32.30850
2       5 29.30842
3      50 29.86611
4     500 29.64647

# 协变量放在前面，分析控制协变量时，自变量的影响
> fit <- aov(weight~gesttime+dose, data = litter)
> summary(fit)
            Df Sum Sq Mean Sq F value  Pr(>F)   
gesttime     1  134.3  134.30   8.049 0.00597 ** # 怀孕时间与出生幼崽的体重有关
dose         3  137.1   45.71   2.739 0.04988 * # 控制怀孕时间，药物剂量与出生体重相关
Residuals   69 1151.3   16.69                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

two-way ANOVA

two-way ANOVA，双因素方差分析，两个自变量

ToothGrowth数据集

自变量：喂食种类supp、剂量dose

响应变量：牙齿长度len

> head(ToothGrowth)
   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

# 将自变量转化为因子类型
> class(ToothGrowth$supp)
[1] "factor"
> class(ToothGrowth$dose)
[1] "numeric"
> ToothGrowth$dose <- factor(ToothGrowth$dose)
> class(ToothGrowth$dose)
[1] "factor"

# 频数统计二联表
> xtabs(~dose+supp, data = ToothGrowth)
     supp
dose  OJ VC
  0.5 10 10
  1   10 10
  2   10 10

# 分组统计均值
> aggregate(ToothGrowth$len, by=list(ToothGrowth$dose, ToothGrowth$supp), FUN=mean)
  Group.1 Group.2     x
1     0.5      OJ 13.23
2       1      OJ 22.70
3       2      OJ 26.06
4     0.5      VC  7.98
5       1      VC 16.77
6       2      VC 26.14

# 二者存在交互，用*表示所有情况，即A+B+A:B
> fit <- aov(len~supp*dose, data = ToothGrowth)
> summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
supp         1  205.4   205.4  15.572 0.000231 *** # 结果显著，supp和dose对len都有影响
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *   # 交互项也显著
Residuals   54  712.1    13.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 可使用HH包的函数将dose和supp的交互效应可视化
> library(HH)
> interaction.plot(ToothGrowth$dose, ToothGrowth$supp, ToothGrowth$len)

多元方差分析

当因变量不止一个时，可以使用多元方差分析

MASS包UScereal数据集

自变量：货架shelf

因变量：谷物所含卡路里calories、脂肪fat、糖sugars

> library(MASS)
> attach(UScereal)

# aggregate()只接受一个值，cbind()连接
> aggregate(calories, fat, sugars, by=list(shelf), FUN=mean)
Error in mean.default(X[[i]], ...) : 'trim'必需是长度必需为一的数值
In addition: Warning message:
In if (na.rm) x <- x[!is.na(x)] :
  the condition has length > 1 and only the first element will be used
> aggregate(cbind(calories, fat, sugars), by=list(shelf), FUN=mean)
  Group.1 calories       fat    sugars
1       1 119.4774 0.6621338  6.295493
2       2 129.8162 1.3413488 12.507670
3       3 180.1466 1.9449071 10.856821

# 自变量转化为因子
> shelf <- factor(shelf)

# +号只能写右边，同样cbind()连接
> fit <- manova(cbind(calories, fat, sugars)~shelf)
> summary(fit)
          Df Pillai approx F num Df den Df    Pr(>F)    
shelf      2 0.4021   5.1167      6    122 0.0001015 ***
Residuals 62                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


# 输出单变量结果
> summary.aov(fit)
 Response calories :
            Df Sum Sq Mean Sq F value    Pr(>F)    
shelf        2  50435 25217.6  7.8623 0.0009054 ***
Residuals   62 198860  3207.4                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response fat :
            Df Sum Sq Mean Sq F value  Pr(>F)  
shelf        2  18.44  9.2199  3.6828 0.03081 *
Residuals   62 155.22  2.5035                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response sugars :
            Df  Sum Sq Mean Sq F value   Pr(>F)   
shelf        2  381.33 190.667  6.5752 0.002572 **
Residuals   62 1797.87  28.998                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Student’s t test

学生 t 检验（student’s t test，也称 t-test）是一种统计检验，用于比较两组的均值。在假设检验中常用于确定某个过程、处理是否对目标真的有显著影响，或者两组是否有显著差异。

比如：想知道鸢尾花的平均花瓣长度是否和花的种类有关，收集两种鸢尾花的花瓣长度数据，每组测量25次，可使用 t 检验检测两组之间的差异。

零假设 H₀：两组均值之间的真实差异为零，即没有差异

备择假设H_a：真实差异不为零

适用情况： t 检验只能用于比较两组的均值（也称成对比较），如果想比较两组以上，用 ANOVA 检验或 post-hoc test。

t 检验假设

数据之间相互独立
数据满足（或近似满足）正态分布
方差同质性：被比较的组的方差相似

若数据不满足要求，可以尝试非参数检验（nonparametric test），如 Wilcoxon Signed-Rank test。

t 检验类型

one-sample，two-sample，paired t-test？

paired t-test，配对样本均值检验：组来自同一类别（如某种处理的前后测量）
two-sample，两独立样本均值检验：组来自两个不同类别（如两种物种，来自两个城市的人），也即 independent t-stet
one-sample，单样本均值检验：一组和标准值相比较（如将液体的 ph 与中性 ph7 相比）

one-tailed，two-tailed t-test？

two-tailed，双侧 t 检验：只关注两组是否有差异
one-tailed，单侧 t 检验：需要知道两组的均值谁大谁小

如想知道鸢尾花的平均花瓣长度是否和花的种类有关，观测来自两种花，故使用 two-sample t-test；只关注是否有差异，故使用 two-tailed t-test。

t 检验方程

$t=\frac{\bar{x_{1}}-\bar{x_{2}}}{\sqrt{s^{2}(\frac{1}{n_{1} }+ \frac{1}{n_{2} })} }$

式中，t 为 t 值，x1，x2是两组的均值，s两组的合并标准误差，n1，n2是每组的观测个数

t 值大，则两组的均值差异大于合并标准误差，表明两组的差异显著

实例

# 检验鸢尾花花瓣长度是否与花的种类有关
# Petal.Length：花瓣长度
# Species：花种类

> flower.data <- read.csv('flower.data.csv', row.names = 1)
> class(flower.data$Species)
[1] "character"

# 转换为因子
> flower.data$Species <- factor(flower.data$Species)
> class(flower.data$Species)
[1] "factor"

> head(flower.data)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.7         3.2          1.3         0.2  setosa
3          5.0         3.6          1.4         0.2  setosa
4          4.6         3.4          1.4         0.3  setosa
5          4.4         2.9          1.4         0.2  setosa
6          5.4         3.7          1.5         0.2  setosa
> table(flower.data$Species)

   setosa virginica 
       25        25 

# two-sample two-tailed t test
> t.test(Petal.Length~Species, flower.data, paired=F, alternative='two.sided')

	Welch Two Sample t-test

data:  Petal.Length by Species
t = -33.719, df = 30.196, p-value < 2.2e-16
alternative hypothesis: true difference in means between group setosa and group virginica is not equal to 0
95 percent confidence interval:
 -4.331287 -3.836713
sample estimates:
   mean in group setosa mean in group virginica 
                  1.456                   5.540

t.test()输出内容解释

data：对比项
t值：-33.719，虽然是负的，但是大多数情况下，只关注其绝对值
df：自由度，30.196，自由度与样本大小有关，表示有多少数据点可用于比较，自由度越大，统计检验越有效
p值：2.2e-16，偶然得到这样的t值的概率
alternative hypothesis：备择假设，setosa 和 virginica 组的真实均值差异不为0
95% confidence interval：95% 置信区间，真实的均值差在 95% 可能性下所处的区间，即两组的真实均值差有 95% 可能性处于[-4.331287, -3.836713]区间内；可以更改区间，但 95% 常用
mean：各组的均值

由以上结果，样本数据的均值差 = 1.456 - 5.540 = -4.084，置信区间显示均值的真实差异在[-4.331287, -3.836713]之间。所以，95% 的情况下均值的真实差异不为0。p值远小于 0.05，所以可以拒绝原假设（均值差异为0），即可以说，均值的真实差异显著不为0。

t 检验总结方式

展示重要的信息：t 值，p 值，自由度

也可以带上均值（mean）和标准差差（standard deviation）

> flower.data %>%
  dplyr::group_by(Species) %>%
  dplyr::summarize(mean_length=mean(Petal.Length),sd_length=sd(Petal.Length))
# A tibble: 2 x 3
  Species   mean_length sd_length
  <fct>           <dbl>     <dbl>
1 setosa           1.46     0.206
2 virginica        5.54     0.569

功效分析

power analysis，可以帮助在给定置信度的情况下，判断检测到给定效应值所需的样本量。反过来，也可以在给定置信度水平情况下，计算某样本量内能检测到给定效应值的概率。

理论基础

样本大小n：实验设计中每种条件/组中观测的数目
显著性水平p（也称alpha）：Type I error probability。也可把他看做是发现效应不发生的概率
功效1-p：1 minus Type II error probability。可把他看成是真实效应发生的概率
效应值ES：在备择假设或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法

功效分析中研究设计的四个基本量，知三求一

R中使用pwr包进行功效分析，根据不同的假设检验选择不同的函数

线性回归的功效分析：pwr.f2.test()

假设要研究老板风格对员工满意度的影响，是否超过薪水和消费对员工满意度的影响。老板风格用4个变量来估计，薪水和小费与3个变量有关。已有研究表明，薪水和小费能够解释约30%员工满意度的方差，而从现实出发，老板风格至少能解释35%的方差，假定显著性水平为0.05，在90%的置信度水平下，需要多少受试者，才能得到这样的方差贡献率？

# Linear Models
# sig.level：显著性水平，默认0.05；power：功效水平；u, v：分子自由度和分母自由度；f2：效应值
> pwr.f2.test(u=3, f2=0.0769, sig.level=0.05, power=0.90)

     Multiple regression power calculation 

              u = 3
              v = 184.2426 # 至少需要185个样本
             f2 = 0.0769   # f2通过计算得来
      sig.level = 0.05
          power = 0.9

方差分析：pwr.anova.test()

假设2组样品做单因素方差分析，要达到0.9的功效，效应值为0.25，并选择0.05的显著性水平，求每一组的样本量

# ANOVA
# k：组数；n：每组的样本大小；f：效应值；sig.level：显著性水平；power：功效
> pwr.anova.test(k=2,f=.25,sig.level=.05,power=.9)

     Balanced one-way analysis of variance power calculation 

              k = 2
              n = 85.03128 # 每组需要86个样本，总样本大小×2
              f = 0.25
      sig.level = 0.05
          power = 0.9

NOTE: n is number in each group

t检验：pwr.t.test()

# t tests
# n：每组样本大小；d：效应值，即标准化的均值之差
# sig.level：显著性水平；power：功效；type：t检验类型（单样本、双样本或相依样本）
# alternative：统计检验为双侧还是单侧检验
> pwr.t.test(d=.8, sig.level=.05,power=.9, type="two.sample",  
           alternative="two.sided")

     Two-sample t test power calculation 

              n = 33.82555
              d = 0.8
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

> pwr.t.test(n=20, d=.5, sig.level=.01, type="two.sample", 
           alternative="two.sided")

     Two-sample t test power calculation 

              n = 20
              d = 0.5
      sig.level = 0.01
          power = 0.1439551
    alternative = two.sided

NOTE: n is number in *each* group

广义线性模型

前面所学的线性回归和方差分析基于正态分布的假设

许多变量并不是正态分布

install.packages('robust')
data(breslow.dat, package='robust')
summary(breslow.dat)

fit <- glm(sumY~Base+Age+Trt, data=breslow.dat, family = poisson(link = 'log'))
summary(fit)
coef(fit)
exp(coef(fit))

R语言作图

条形图

# 读数据
x <- read.csv('../Rdata/homo_length.csv') # 可加参数 header = FALSE，不读表头
# 取所需数据
xx
xx <- x[1:24,]
xx$chr # $查看
xx$length
# 画图
barplot(height = xx$length) # 先画基本图，避免报错
barplot(height = xx$length, names.arg = xx$chr) # 每个条形加标签
barplot(height = xx$length, names.arg = xx$chr, las = 2) # 改变标尺数字方向
yanse <- sample(colours(),24,replace = F) # 设置随机颜色
barplot(height = xx$length, names.arg = xx$chr, las = 2, col = yanse)
barplot(height = xx$length, names.arg = xx$chr, las = 2, col = yanse, border = F) # 取消条形图描边
barplot(height = xx$length, names.arg = xx$chr, las = 3, col = yanse, border = F, main = 'Human Chromosome Length Distribution', xlab = 'Chromosome', ylab = 'Length') # 设置图表标题、坐标轴标题

分组条形图

# 读数据
x <- read.csv(file = '../Rdata/sv_distrubution.csv', row.names = 1) # 第一列作为行名
tail(x) # 查看最末或最前几行
head(x)
# 画图
barplot(as.matrix(x)) # x为数据框，转换为矩阵
barplot(t(as.matrix(x))) # 转置一下
barplot(t(as.matrix(x)), col = rainbow(4), beside = T) # 默认堆砌，可以转换为分组
colnames(x) # 矩阵每列名称，即表头
barplot(t(as.matrix(x)), col = rainbow(4), legend.text = colnames(x)) #将图例设置为列名
barplot(t(as.matrix(x)), col = rainbow(4), legend.text = colnames(x), ylim = c(0,800)) # 设置y轴范围
barplot(t(as.matrix(x)), col = rainbow(4), legend.text = colnames(x), ylim = c(0,800), main = 'SV Distribution', ylab = 'Numbers', xlab = 'Chromosomes')

直方图

# 读数据，跳过前7行，字符串取消用''显示
x <- read.table('../Rdata/H37Rv.gff', sep = '\t', header = F, skip = 7, quote = '')
View(x)
head(x$V3=='gene') #取类型为gene的行
[1] FALSE  TRUE FALSE  TRUE FALSE  TRUE
x[x$V3=='gene',]
x <- x[x$V3=='gene',]
x <- abs(x$V5-x$V4+1) # 计算基因长度，保险起见加绝对值
length(x)
[1] 3978
range(x) # 范围
[1]    61 12456
hist(x)
hist(x, breaks = 80) # 分成80个
hist(x, breaks = c(0,500,1000,1500,2000,2500,15000)) # 也可以手动设置范围
hist(x, breaks = 80, freq = F) # 默认以频数显示，可以改成频率
hist(x, breaks = 80, freq = F, density = T)
hist(x, breaks = 80, freq = F, density = T, col = 'pink')
hist(x, breaks = 80, freq = F, col = 'pink', border = F, main = 'histogram of Gene Length', xlab = 'Gene Length (bp)')

h <- hist(x, breaks = 80, freq = F, col = 'pink', border = F, main = 'histogram of Gene Length', xlab = 'Gene Length (bp)')
class(h) # h包含了此直方图的所有信息
[1] "histogram"
h$xname
[1] "x"
h$breaks
 [1]     0   200   400   600   800  1000  1200  1400  1600  1800  2000  2200  2400  2600  2800  3000
[17]  3200  3400  3600  3800  4000  4200  4400  4600  4800  5000  5200  5400  5600  5800  6000  6200
[33]  6400  6600  6800  7000  7200  7400  7600  7800  8000  8200  8400  8600  8800  9000  9200  9400
[49]  9600  9800 10000 10200 10400 10600 10800 11000 11200 11400 11600 11800 12000 12200 12400 12600
h$counts
 [1] 120 548 578 601 571 430 356 268 165  84  55  43  27  20  24  12  12  13   7   5   3   5   4   4
[25]   5   1   1   1   2   0   0   4   1   1   0   0   0   2   0   0   0   0   0   0   0   0   1   1
[49]   0   1   0   0   0   0   0   1   0   0   0   0   0   0   1
h$density
 [1] 1.508296e-04 6.887883e-04 7.264957e-04 7.554047e-04 7.176973e-04 5.404726e-04 4.474610e-04
 [8] 3.368527e-04 2.073906e-04 1.055807e-04 6.913022e-05 5.404726e-05 3.393665e-05 2.513826e-05
[15] 3.016591e-05 1.508296e-05 1.508296e-05 1.633987e-05 8.798391e-06 6.284565e-06 3.770739e-06
[22] 6.284565e-06 5.027652e-06 5.027652e-06 6.284565e-06 1.256913e-06 1.256913e-06 1.256913e-06
[29] 2.513826e-06 0.000000e+00 0.000000e+00 5.027652e-06 1.256913e-06 1.256913e-06 0.000000e+00
[36] 0.000000e+00 0.000000e+00 2.513826e-06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[43] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.256913e-06 1.256913e-06 0.000000e+00
[50] 1.256913e-06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.256913e-06
[57] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.256913e-06

散点图

# 读数据
m <- read.table('../Rdata/prok_representative.txt', sep = '\t')
View(m)
gene_number <- m$V4
gnome_size <- m$V2

# 设置散点类型(0-25)，散点大小
plot(gnome_size, gene_number,pch=16,cex=0.8,xlab = 'Gnome Size',ylab = 'Genes')

# 线性回归
fit <- lm(gene_number ~ gnome_size)
fit

Call:
lm(formula = gene_number ~ gnome_size)

Coefficients:
(Intercept)   gnome_size  
      286.6        843.7  

# 画线
abline(fit, col='red')
# 可以直接根据回归结果添加文本
text(12,6000,labels='y=843.691x+286.598\nR2=0.97')

# 也可以设置变量取到相应值，这样不用每次都手动敲
c # R2存储位置
[1] 0.9676204
round(summary(fit)$adj.r.squared, 2)
[1] 0.97
rr <- round(summary(fit)$adj.r.squared, 2)
intercept <- round(summary(fit)$coefficients[1], 2) # 截距存储位置
intercept
[1] 286.6
slope <- round(summary(fit)$coefficients[2], 2) # 斜率存储位置
slope
[1] 843.69
eq <- bquote( atop( "y = " * .(slope) * " x + " * .(intercept), R^2 == .(rr) ) )
text(12,6000,eq)

饼图

# 读文件
read.table('../Rdata/Species.txt')
        V1             V2       V3
1  Viruses         noname 31.42034
2 Bacteria  Bacteroidetes 31.27472
3 Bacteria  Fibrobacteres 24.74666
4 Bacteria     Firmicutes  8.06092
5 Bacteria   Spirochaetes  2.32585
6  Archaea  Euryarchaeota  1.42253
7 Bacteria Proteobacteria  0.74899
x <- read.table('../Rdata/Species.txt')

# 画饼图
pie(x$V3)
# 设置数据标签
lab <- paste(x$V2,'\n',x$V3,'%')
lab
[1] "noname \n 31.42034 %"        "Bacteroidetes \n 31.27472 %" "Fibrobacteres \n 24.74666 %"
[4] "Firmicutes \n 8.06092 %"     "Spirochaetes \n 2.32585 %"   "Euryarchaeota \n 1.42253 %" 
[7] "Proteobacteria \n 0.74899 %"
# 设置颜色、饼图半径、标签大小
# pie数据个数多的话，标签容易堆在一起，可以调小字体，还是重叠可以用plotrix的pie.labels()函数
pie(x$V3, labels = lab, col = rainbow(nrow(x)))
pie(x$V3, labels = lab, col = rainbow(nrow(x)), radius = 1, cex = 0.8)
# 3D饼图需要plotrix库
install.packages("plotrix")
library(plotrix)s
# pie换成pie3D即可，可以设置explode分割距离
pie3D(x$V3, labels = lab, col = rainbow(nrow(x)), radius = 1, cex = 0.8)
pie3D(x$V3, labels = lab, col = rainbow(nrow(x)), radius = 1, cex = 0.8, explode = 0.2)

# pie3D.labels()解决标签堆叠问题
pieplot3d <- pie(x$V3, col = rainbow(nrow(x)), radius = 1, explode = 0.1)
pie3D.labels(pieplot3d,labels = lab,labelcex = 0.8,height = 0.3,labelrad = 1.4)

# pie.labels()解决标签堆叠问题
pieplot3d <- pie(x$V3, labels = '',col = rainbow(nrow(x)), radius = 1)
ang <- floating.pie(0,0,x$V3) # 标签分散的角度，可以设置成数据的百分比
pie.labels(x=0,y=0, angles = ang, labels=lab, radius=1.2, minangle = 0.17) # 最小角度设置为10/180*pie

箱线图

箱线图展示了数据的5个指标（最小值、下四分位数、中位数、上四分位数、最大值）和离群点。

boxplot(mtcars$mpg)
fivenum(mtcars$mpg) # 给出数据的最小值、下四分位数、中位数、上四分位数、最大值
[1] 10.40 15.35 19.20 22.80 33.90

# x ~ y，y为分组因子，多个y用:隔开
boxplot(mpg ~ cyl, data = mtcars)
boxplot(hp ~ cyl, data = mtcars)
boxplot(len ~ dose, data = ToothGrowth)
boxplot(len ~ dose:supp, data = ToothGrowth)
boxplot(mpg ~ cyl, data = mtcars, col = c('pink', 'orange', 'red')) # 颜色
boxplot(mpg ~ cyl, data = mtcars, col = c('pink', 'orange', 'red'), pch = 16, cex = 0.8) # 离散点形状、大小
boxplot(mpg ~ cyl, data = mtcars, col = c('pink', 'orange', 'red'), pch = 16, cex = 0.8, xlab = 'number of cylinders', ylab = 'miles/gallon')

boxplot(len ~ dose, data = ToothGrowth, notch = T) # 显示缺口
boxplot(len ~ dose, data = ToothGrowth, width = c(0.5,1,2)) # 指定箱子宽度
boxplot(len ~ dose, data = ToothGrowth, varwidth = T) # 箱子宽度与组中观察数的平方根成正比
boxplot(len ~ dose, data = ToothGrowth, boxwex = 0.5) # 所有箱子宽度等比缩放
boxplot(len ~ dose, data = ToothGrowth, staplewex = 0.5) # 最大最小值横线的宽度
boxplot(len ~ dose:supp, data = ToothGrowth, sep = ':') # 多个分组因子y的分隔符号
boxplot(len ~ dose:supp, data = ToothGrowth, sep = ':', lex.order = T) # 箱子排序

# 读数据，第一列作为行名
x <- read.csv('../Rdata/gene_expression.csv', header = T, row.names = 1)
head(x)
         A1    A2    A3    B1    B2    B3
gene1  0.33  0.13  0.30  0.48  0.37  1.73
gene2  0.00  0.00  0.30  0.54  0.43  0.35
gene3 10.27 10.05  5.51 30.95 15.85 25.72
gene4  3.74  4.44  4.24  0.24  0.06  1.04
gene5  9.66 10.05 15.57  6.83  2.66  5.89
gene6  0.50  0.64  0.18  0.72  1.24  2.60
boxplot(x) # 可以直接绘制
boxplot(x, outline = F) # 去除离群点

# 也可以用reshape2函数重新组合数据，用x ~ y的方式，推荐这种
library(reshape2)
x <- melt(x)
No id variables; using all as measure variables
head(x)
  variable value
1       A1  0.33
2       A1  0.00
3       A1 10.27
4       A1  3.74
5       A1  9.66
6       A1  0.50
boxplot(value ~ variable, data = x, outline = F)

热图

生信分析中常见的图，根据颜色的深浅反应数值的大小，可用于基因表达的可视化和聚类。

绘制热图的对象是数值矩阵。

绘制热图的方法：

heatmap()
gplot包heatmap.2()
pheatmap包pheatmap()
ComplexHeatmap包

# 不是数值矩阵要先转换
class(mtcars)
[1] "data.frame"
heatmap(mtcars)
Error in heatmap(mtcars) : 'x' must be a numeric matrix
heatmap(as.matrix(mtcars))

# 如果数据差距过大，热图无法分辨，可以加入scale参数按row或col进行标准化
heatmap(as.matrix(mtcars), scale = 'col')

# 读数据
m <- read.csv('../Rdata/heatmap.csv', row.names = 1)
x <- as.matrix(m)
heatmap(x)
heatmap(t(x)) # 可转置互换坐标轴

# 添加颜色，但是无渐变不好看
heatmap(x, col=c('green','blue'))

# 生成渐变色
colorRampPalette(c('red', 'black', 'green'))(nrow(x))
 [1] "#FF0000" "#ED0000" "#DB0000" "#CA0000" "#B80000" "#A70000" "#950000" "#830000" "#720000"
[10] "#600000" "#4F0000" "#3D0000" "#2B0000" "#1A0000" "#080000" "#000800" "#001A00" "#002B00"
[19] "#003D00" "#004F00" "#006000" "#007200" "#008300" "#009500" "#00A700" "#00B800" "#00CA00"
[28] "#00DB00" "#00ED00" "#00FF00"
yanse <- colorRampPalette(c('red', 'black', 'green'))(nrow(x))
heatmap(x, col = yanse)
# 添加边色
heatmap(x, col = yanse, ColSideColors = colorRampPalette(c('red', 'black', 'green'))(ncol(x)), RowSideColors = yanse)

# 取消按列或按行聚类
heatmap(x, col = yanse, Colv = NA, Rowv = NA)

# heatmap.2()绘制热图
library(gplots)
m <- as.matrix(read.csv('../Rdata/heatmap.csv', row.names = 1))
heatmap.2(m)
heatmap.2(m, key = F) # 是否显示图例
heatmap.2(m, symkey = T) # 图例是否对称
heatmap.2(m, symkey = F)
heatmap.2(m, symkey = F, density.info = 'none')
heatmap.2(m, symkey = F, trace = 'none') # 是否显示图例
heatmap.2(m, symkey = F, tracecol = 'green') # 竖线颜色
heatmap.2(m, dendrogram = 'none') # 是否聚类
heatmap.2(m, dendrogram = 'row')
heatmap.2(m, dendrogram = 'col')

# pheatmap()
library(pheatmap)
pheatmap(m)
pheatmap(m, border_color = 'white')
# 加入聚类类型
rep(c('N1','T1','N2','T2'), each=5)
 [1] "N1" "N1" "N1" "N1" "N1" "T1" "T1" "T1" "T1" "T1" "N2" "N2" "N2" "N2" "N2" "T2" "T2" "T2"
[19] "T2" "T2"
factor(rep(c('N1','T1','N2','T2'), each=5))
 [1] N1 N1 N1 N1 N1 T1 T1 T1 T1 T1 N2 N2 N2 N2 N2 T2 T2 T2 T2 T2
Levels: N1 N2 T1 T2
annotation_col <- data.frame(factor(rep(c('N1','T1','N2','T2'), each=5)))
row.names(annotation_col) <- colnames(m)
colnames(annotation_col) <- 'Cell Type'
head(annotation_col)
      Cell Type
N.GD1        N1
N.GD2        N1
N.GD3        N1
N.GD4        N1
N.GD5        N1
T.GD1        T1

图表设置

图表参数

par函数，Parameter简写，设置绘图中的参数。

head(par())
$xlog
[1] FALSE

$ylog
[1] FALSE

$adj
[1] 0.5

$ann
[1] TRUE

$ask
[1] FALSE

$bg
[1] "white"

?par() # 查看帮助文档
# par()中自带的为默认值，可以更改，如
par(pch=16) # 更改点样式，之后绘图均按此值，重启R后恢复到默认值
# 为防止默认值更改，存储可更改的参数
opar <- par(no.readonly = T)
par(opar) # 恢复默认值

# 具体设置查看帮助文档
plot(women$height,women$weight,type = "b",pch=16,lty=2,lwd=1.2,las=1,
     main = "This is main title", sub = "This is sub title",xlab = "x label",ylab = "y label",
     xlim = c(50,100),ylim = c(100,200),
     adj=0.5,bty="]",tck=-0.03,
     col="red",fg="blue",col.main="green",col.sub='pink',col.lab='orange',col.axis='purple',
     cex=1.5,cex.lab=1.2,cex.axis=1.2,cex.main=1.5)

图表颜色

颜色选择标准：

对比明显
低调内敛
色盲友好

colors()
  [1] "white"                "aliceblue"            "antiquewhite"        
  [4] "antiquewhite1"        "antiquewhite2"        "antiquewhite3"       
  [7] "antiquewhite4"        "aquamarine"           "aquamarine1" 
...

# 随机抽取颜色
sample(colors(),5)
[1] "cadetblue4"      "palegreen4"      "hotpink2"        "lightsteelblue2" "gray12" 

rgb(0.1, 0.2, 0.3)
[1] "#1A334D"
hsv(0.1, 0.2, 0.3) # 色调（H），饱和度（S），明度（V）
[1] "#4D463D"
hcl(0.1, 0.2, 0.3) # 色调（H），色度（S），亮度（L）
[1] "#020101"
rainbow(5) # 彩虹色
[1] "#FF0000" "#CCFF00" "#00FF66" "#0066FF" "#CC00FF"

heat.colors(5)
[1] "#FF0000" "#FF5500" "#FFAA00" "#FFFF00" "#FFFF80"
terrain.colors(5)
[1] "#00A600" "#E6E600" "#EAB64E" "#EEB99F" "#F2F2F2"
topo.colors(5)
[1] "#4C00FF" "#004CFF" "#00E5FF" "#00FF4D" "#FFFF00"
cm.colors(5)
[1] "#80FFFF" "#BFFFFF" "#FFFFFF" "#FFBFFF" "#FF80FF"
grey.colors(5)
[1] "#4D4D4D" "#888888" "#AEAEAE" "#CCCCCC" "#E6E6E6"

# R包RColorBrewer带了很多颜色
library(RColorBrewer)
display.brewer.all() # 显示所有颜色集
display.brewer.pal(name = 'Set3', n = 5) # 显示所选颜色
brewer.pal(name = 'Set3', n = 5) # 设置颜色
[1] "#8DD3C7" "#FFFFB3" "#BEBADA" "#FB8072" "#80B1D3"
pie(1:5, col = brewer.pal('Set3', n = 5)) # 使用时传递给col参数即可

韦恩图

一般最少2组，最多5组，再多图很复杂没有意义。

韦恩图画法：

gplots包venn函数（简单，无颜色）
venneuler包中的venneuler函数（安装失败）
venn包中的venn函数（简单，有颜色，推荐）
vennDiagram包中的venn.Diagram函数（有颜色，图没有上个好看，必须输出文件，无法在plots窗口预览）

# 读数据
listA <- read.csv('../Rdata/genes_list_A.txt', header = F)
listB <- read.csv('../Rdata/genes_list_B.txt', header = F)
listC <- read.csv('../Rdata/genes_list_C.txt', header = F)
listD <- read.csv('../Rdata/genes_list_D.txt', header = F)
listE <- read.csv('../Rdata/genes_list_E.txt', header = F)
A <- listA$V1
B <- listB$V1
C <- listC$V1
D <- listD$V1
E <- listE$V1
length(A);length(B);length(C);length(D);length(E)
[1] 562
[1] 474
[1] 2273
[1] 1735
[1] 594

# gplots包
library(gplots)
# 绘图简单，没有颜色
venn(list(A,B,C,D,E))

# vennDiagram包
library(VennDiagram)
# 需指明输出文件，无法预览
venn.diagram(list(sampleC=C, sampleD=D), filename = 'venn.png', fill = c('yellow', 'cyan'), cex = 1.5)
[1] 1
# 不输出log文件
venn.diagram(list(sampleC=C, sampleD=D), filename = 'venn2.png', fill = c('yellow', 'cyan'), cex = 1.2, disable.logging = T)
INFO [2022-02-01 17:06:37] [[1]]
INFO [2022-02-01 17:06:37] list(sampleC = C, sampleD = D)
INFO [2022-02-01 17:06:37] 
INFO [2022-02-01 17:06:37] $filename
INFO [2022-02-01 17:06:37] [1] "venn2.png"
INFO [2022-02-01 17:06:37] 
INFO [2022-02-01 17:06:37] $fill
INFO [2022-02-01 17:06:37] c("yellow", "cyan")
INFO [2022-02-01 17:06:37] 
INFO [2022-02-01 17:06:37] $cex
INFO [2022-02-01 17:06:37] [1] 1.2
INFO [2022-02-01 17:06:37] 
INFO [2022-02-01 17:06:37] $disable.logging
INFO [2022-02-01 17:06:37] T
INFO [2022-02-01 17:06:37] 
[1] 1
# 三维
venn.diagram(list(sampleC=C, sampleD=D, sampleA=A), filename = 'venn3.png', fill = c('yellow', 'cyan', 'blue'), cex = 1.2, disable.logging = T)
# 四维
venn.diagram(list(sampleC=C, sampleD=D, sampleA=A, sampleB=B), filename = 'venn4.png', fill = c('yellow', 'cyan', 'blue', 'red'), cex = 1.2, disable.logging = T)
# 五维
venn.diagram(list(sampleC=C, sampleD=D, sampleA=A, sampleB=B, sampleE=E), filename = 'venn5.png', fill = c('yellow', 'cyan', 'blue', 'red', 'pink'), cex = 1.2, disable.logging = T)

# venn包
venn(list(A,B,C,D,E))
# 不显示方框
venn(list(A,B,C,D,E), box = F)
# 改变边线颜色
venn(list(A,B,C,D,E), box = F, col='red')
# 加标签
venn(list(A,B,C,D,E), box = F, snames = 'Sample A,Sample B,Sample C,Sample D,Sample E')
# 加颜色
venn(list(A,B,C,D,E), box = F, snames = 'Sample A,Sample B,Sample C,Sample D,Sample E', zcolor = 'red, green, purple, blue, cyan')
# 如果可以，强制形状为椭圆
venn(list(A,B,C,D,E), box = F, snames = 'Sample A,Sample B,Sample C,Sample D,Sample E', zcolor = 'red, green, purple, blue, cyan', ellipse = T)

曼哈顿图

曼哈顿图是用二维散点图展示大量数据的方法，主要用于全基因组关联分析。

library(qqman)
library(RColorBrewer)

# 此包内置数据集，导入数据时也要是这样的格式
gwasResults
      SNP CHR  BP            P
1     rs1   1   1 0.9148060435
2     rs2   1   2 0.9370754133
3     rs3   1   3 0.2861395348
...

manhattan(gwasResults)

# 加标题，设置y轴范围，散点大小，轴字体大小
manhattan(gwasResults, main = 'Manhattan Plot', ylim = c(0,6), cex = 0.6, cex.axis = 0.9)

# 控制两根限制条：-log10(1e-5)和-log10(1e-8)
manhattan(gwasResults, main = 'Manhattan Plot', ylim = c(0,6), cex = 0.6, suggestiveline = FALSE, genomewideline = FALSE)

# 加颜色
manhattan(gwasResults, main = 'Manhattan Plot', ylim = c(0,6), cex = 0.6, cex.axis=0.9, col = c('blue4', 'orange3'))
yanse <- brewer.pal('Set3', n=4)
manhattan(gwasResults, main = 'Manhattan Plot', ylim = c(0,6), cex = 0.6, cex.axis=0.9, col = yanse, chrlabs = c(1:20, 'P', 'Q'))

# 染色体标签
manhattan(gwasResults, main = 'Manhattan Plot', ylim = c(0,6), cex = 0.6, cex.axis=0.9, col = c('blue4', 'orange3'), chrlabs = c(1:20, 'P', 'Q'))

# 局部显示显著的染色体
manhattan(subset(gwasResults,CHR==3))

# 输出感兴趣的snps
snpsOfInterest
  [1] "rs3001" "rs3002" "rs3003" "rs3004" "rs3005" "rs3006" "rs3007" "rs3008" "rs3009" "rs3010"
 [11] "rs3011" "rs3012" "rs3013" "rs3014" "rs3015" "rs3016" "rs3017" "rs3018" "rs3019" "rs3020"
 [21] "rs3021" "rs3022" "rs3023" "rs3024" "rs3025" "rs3026" "rs3027" "rs3028" "rs3029" "rs3030"
 [31] "rs3031" "rs3032" "rs3033" "rs3034" "rs3035" "rs3036" "rs3037" "rs3038" "rs3039" "rs3040"
 [41] "rs3041" "rs3042" "rs3043" "rs3044" "rs3045" "rs3046" "rs3047" "rs3048" "rs3049" "rs3050"
 [51] "rs3051" "rs3052" "rs3053" "rs3054" "rs3055" "rs3056" "rs3057" "rs3058" "rs3059" "rs3060"
 [61] "rs3061" "rs3062" "rs3063" "rs3064" "rs3065" "rs3066" "rs3067" "rs3068" "rs3069" "rs3070"
 [71] "rs3071" "rs3072" "rs3073" "rs3074" "rs3075" "rs3076" "rs3077" "rs3078" "rs3079" "rs3080"
 [81] "rs3081" "rs3082" "rs3083" "rs3084" "rs3085" "rs3086" "rs3087" "rs3088" "rs3089" "rs3090"
 [91] "rs3091" "rs3092" "rs3093" "rs3094" "rs3095" "rs3096" "rs3097" "rs3098" "rs3099" "rs3100"

# 高亮显示感兴趣的snps
manhattan(gwasResults, highlight = snpsOfInterest)

# 标注满足p值的rs
manhattan(gwasResults, annotatePval = 0.001)
# 默认只显示最顶部的rs，可全部显示
manhattan(gwasResults, annotatePval = 0.001, annotateTop = FALSE)

火山图

散点图的一种，主要用于基因表达的分析中。

自带散点图

# 读数据
m <- read.csv('../Rdata/Deseq2.csv', header = T, row.names = 1)
head(m)
                    baseMean log2FoldChange      lfcSE       stat       pvalue         padj
ENSG00000000003  708.6021697    -0.38125397 0.10065597 -3.7876937 1.520521e-04 0.0012815112
ENSG00000000419  520.2979006     0.20681259 0.11222180  1.8428915 6.534485e-02 0.1962308610
ENSG00000000457  237.1630368     0.03792034 0.14345322  0.2643394 7.915184e-01 0.9112208706
ENSG00000000460   57.9326331    -0.08816367 0.28716771 -0.3070111 7.588349e-01 0.8946714454
ENSG00000000938    0.3180984    -1.37822703 3.49987280 -0.3937935 6.937335e-01           NA
ENSG00000000971 5817.3528677     0.42640216 0.08831006  4.8284666 1.375884e-06 0.0000181641

# 去除NA
m <- na.omit(m)

# 画图
plot(m$log2FoldChange,m$padj)
plot(m$log2FoldChange,-log10(m$padj), pch=16, cex=0.6, xlim = c(-10,10), ylim = c(0,100))

# 一次性计算-log10(m$padj)并赋值
transform(m,padj=-log10(m$padj))
                   baseMean log2FoldChange      lfcSE         stat       pvalue         padj
ENSG00000000003  708.602170  -0.3812539747 0.10065597 -3.787693677 1.520521e-04 2.892278e+00
ENSG00000000419  520.297901   0.2068125935 0.11222180  1.842891491 6.534485e-02 7.072327e-01
ENSG00000000457  237.163037   0.0379203397 0.14345322  0.264339418 7.915184e-01 4.037634e-02
...

m <- transform(m,padj=-log10(m$padj))

# 将数据分开，为明显上调、明显下调、不显著
down <- m[m$log2FoldChange <= -1,]
up <- m[m$log2FoldChange >= 1,]
mid <- m[m$log2FoldChange -1 & m$log2FoldChange < 1,]

# 分别画图
plot(mid$log2FoldChange, mid$padj, xlim = c(-10,10), ylim = c(0,100), pch = 16, cex = 0.6, col = 'grey', main = 'Gene Expression', xlab = 'log2FodeChange', ylab = '-log10(Qvalue)')

# 低级绘图命令plots加上另外两个
points(up$log2FoldChange,up$padj,pch=16,cex=0.6,col='red')
points(down$log2FoldChange,down$padj,pch=16,cex=0.6,col='green')

ggplot包

# 读数据并处理
m <- read.csv('../Rdata/Deseq2.csv', header = T, row.names = 1)
m <- na.omit(m)
transform(m,padj=-log10(m$padj))
m <- transform(m,padj=-log10(m$padj))

ggplot2作图

逻辑清晰，绘图方便，简介好看，有很多基于ggplot2的绘图包。

图形语法：grammar of graphis，类似于图层的概念。

1.数据（data）

一般为数据框

1	ggplot(data = mtcars)

2.映射（mapping）

可以将size，shape，col等映射为数据框中的值

ggplot(data = mtcars, mapping = aes(x=wt, y=mpg))
ggplot(data = mtcars, mapping = aes(x=wt, y=mpg, size=mpg))
ggplot(data = mtcars, mapping = aes(x=wt, y=mpg, shape=as.factor(cyl)))
ggplot(data = mtcars, mapping = aes(x=wt, y=mpg, col=mpg))

3.几何对象

1	ggplot(data = mtcars, mapping = aes(x=wt, y=mpg)) + geom_point()