03. R <- 基本统计分析Ⅰ - 简要分析、频数表和列联表、基本数学和统计函数

这篇文章主要讲解如何对一些数据进行基本的统计学分析

一、一键式描述性统计分析

针对连续型数据，一般在分析的初始阶段需要获得它的中心趋势、变化性和分布形状。

在这个部分，使用的示例数据为 mtcars

image-20200725175446590

总体计算样本描述性统计量

1. summary & fivenum 一步获得常见统计量

对于一个数据集，可以使用 summary() 来计算常用的描述性统计量：

summary() 提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。

例：

1
2
3
4
# 指定在mtcars中要分析的列
myvars <- c("mpg", "hp", "wt")
# 进行初步分析
summary(mtcars[myvars])

结果：

1
2
3
4
5
6
7
8
> summary(mtcars[myvars])
      mpg              hp              wt       
 Min.   :10.40   Min.   : 52.0   Min.   :1.513  
 1st Qu.:15.43   1st Qu.: 96.5   1st Qu.:2.581  
 Median :19.20   Median :123.0   Median :3.325  
 Mean   :20.09   Mean   :146.7   Mean   :3.217  
 3rd Qu.:22.80   3rd Qu.:180.0   3rd Qu.:3.610  
 Max.   :33.90   Max.   :335.0   Max.   :5.424

fivenum() 可返回图基五数总括（Tukey’s five-number summary，即最小值、下四分位数、中位数、上四分位数和最大值）。

例：

1
2
> fivenum(mtcars$mpg)
[1] 10.40 15.35 19.20 22.80 33.90

2. sapply & apply 通过（自定义）函数进行分析

可以使用apply()函数或sapply()函数计算所选择的任意描述性统计量。对于sapply()函数，其使用格式为：

sapply(x, FUN, options)，x是你的数据框（或矩阵，FUN为一个任意的函数。如果指定了options它们将被传递

给FUN。

  • 可以在这里插入的典型函数有mean()、sd()、var()、min()、max()、median(), length()、range()、quantile()
  • 也可以使用自定义的函数

相关参数：

参数 描述
na.omit TRUE为忽略缺省值

例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 使用预定义函数
> sapply(mtcars[myvars], min)
   mpg     hp     wt 
10.400 52.000  1.513 

# 使用自定义的函数
>  mystats <- 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=n, mean=m, stdev=s, skew=skew, kurtosis=kurt)) 
 }
> myvars <- c("mpg", "hp", "wt") 
> sapply(mtcars[myvars], mystats)
               mpg          hp          wt
n        32.000000  32.0000000 32.00000000
mean     20.090625 146.6875000  3.21725000
stdev     6.026948  68.5628685  0.97845744
skew      0.610655   0.7260237  0.42314646
kurtosis -0.372766  -0.1355511 -0.02271075

3. 第三方包提供的其他方法

若干用户贡献包都提供了计算描述性统计量的函数，其中包括Hmisc、pastecs和psych。由于这些包并未包括在基础安装中，所以需要在首次使用之前先进行安装

Hmisc

Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。

1
2
3
library(Hmisc)
myvars <- c("mpg", "hp", "wt")
describe(mtcars[myvars])

结果：

image-20200725182442001

pastecs

pastecs包中有一个名为stat.desc()的函数，它可以计算种类繁多的描述性统计量。使用格式为：

stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)

参数 描述
x 一个数据框或时间序列
basic 若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。
desc 若desc=TRUE（同样也是默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。
norm 若norm=TRUE（不是默认的），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro-Wilk正态检验结果。
p 使用p值来计算平均数的置信区间（默认置信度为0.95）

例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> library(pastecs) 
> myvars <- c("mpg", "hp", "wt") 
> stat.desc(mtcars[myvars])
                     mpg           hp          wt
nbr.val       32.0000000   32.0000000  32.0000000
nbr.null       0.0000000    0.0000000   0.0000000
nbr.na         0.0000000    0.0000000   0.0000000
min           10.4000000   52.0000000   1.5130000
max           33.9000000  335.0000000   5.4240000
range         23.5000000  283.0000000   3.9110000
sum          642.9000000 4694.0000000 102.9520000
median        19.2000000  123.0000000   3.3250000
mean          20.0906250  146.6875000   3.2172500
SE.mean        1.0654240   12.1203173   0.1729685
CI.mean.0.95   2.1729465   24.7195501   0.3527715
var           36.3241028 4700.8669355   0.9573790
std.dev        6.0269481   68.5628685   0.9784574
coef.var       0.2999881    0.4674077   0.3041285

psych

psych包也拥有一个名为describe()的函数，它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误。

例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> library(psych) 
# 相关先前引入的同名函数被覆盖
Attaching package: 'psych' 
The following object is masked from ‘package:Hmisc’:

    describe

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha
> myvars <- c("mpg", "hp", "wt")
> describe(mtcars[myvars])
    vars  n   mean    sd median trimmed   mad   min    max  range skew kurtosis    se
mpg    1 32  20.09  6.03  19.20   19.70  5.41 10.40  33.90  23.50 0.61    -0.37  1.07
hp     2 32 146.69 68.56 123.00  141.19 77.10 52.00 335.00 283.00 0.73    -0.14 12.12
wt     3 32   3.22  0.98   3.33    3.15  0.77  1.51   5.42   3.91 0.42    -0.02  0.17
  • 最后载入的程序包优先。在这里，psych在Hmisc之后被载入，然后显示了一条信息，提示Hmisc包中的describe()函数被psych包中的同名函数所屏蔽（masked）。键入describe()后，R在搜索这个函数时将首先找 到psych包中的函数并执行它。如果你想改而使用Hmisc包中的版本，可以键入Hmisc::describe(mt)。这个函数仍然在那里。你只是需要给予R更多信息以找到它。

分组计算描述性统计量

在比较多组个体或观测时，关注的焦点经常是各组的描述性统计信息，而不是样本整体的描述性统计信息。同样地，在R中完成这个任务有若干种方法。我们将以获取变速箱类型各水平的描述性统计量开始。

1. 整合(Aggregate)和重塑(reshape)数据

R中提供了许多用来整合（aggregate）和重塑（reshape）数据的强大方法。在整合数据时，往往将多组观测替换为根据这些观测计算的描述性统计量。在重塑数据时，则会通过修改数据的结构（行和列）来决定数据的组织方式。

Transposition 转置

转置（反转行和列）也许是重塑数据集的众多方法中最简单的一个了。使用函数t()即可对一个矩阵或数据框进行转置。对于后者，行名将成为变量（列）名。

例：

1
2
3
4
5
6
7
8
9
10
11
12
> cars <- mtcars[1:5,1:2]
> cars
                   mpg cyl
Mazda RX4         21.0   6
Mazda RX4 Wag     21.0   6
Datsun 710        22.8   4
Hornet 4 Drive    21.4   6
Hornet Sportabout 18.7   8
> t(cars)
    Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive Hornet Sportabout
mpg        21            21       22.8           21.4              18.7
cyl         6             6        4.0            6.0               8.0

使用 aggregate & by 整合

aggregate

aggregate 函数是用来对数据进行分组（折叠）的函数。

在R中使用一个或多个by变量和一个预先定义好的函数来折叠（collapse）数据是比较容易的。调用格式为：

aggregate(x, by, FUN)

参数 描述
x 待折叠的数据对象
by 分组因子。为一个变量名组成的列表，这些变量将被去掉以形成新的观测。R将根据这个列表内的值来进行分组。
FUN 用来计算描述性统计量的标量函数，它将被用来计算新分出组的值。
  • 在使用aggregate()函数的时候，by中的变量必须在一个列表中（即使只有一个变量）。你可以在列表中为各组声明自定义的名称，例如by=list(Group.age=age, Group.gender=gender)。

例：

构建示例用数据集：

1
2
3
4
5
6
7
8
9
10
11
12
13
> name <- c('a','h','b','c','d','e')
> age <- c(18,16,12,14,16,14)
> gender <- c('M','M',"F",'F',"F",'F')
> height <- c(180,170,120,140,160,150)
> frame1 <- data.frame(name=name,age=age,gender=gender,height=height)
> frame1
  name age gender height
1    a  18      M    180
2    b  16      M    170
3    c  12      F    120
4    d  14      F    140
5    e  16      F    160
6    f  14      F    150

单因子分组：在这个例子当中，按照年龄分组，并且对同组求均值。其中无法被mean函数计算的值，将表示为NA

1
2
3
4
5
6
> aggregate(frame1, by = list(Group.age=frame1$age), FUN = mean)
  Group.age name age gender height
1        12   NA  12     NA    120
2        14   NA  14     NA    145
3        16   NA  16     NA    165
4        18   NA  18     NA    180

多因子分组：在这个例子当中，按照年龄和性别分组，并求每组的最大值。

1
2
3
4
5
6
7
> aggregate(frame1, by = list(Group.age=frame1$age, Group.gender=frame1$gender), FUN = max)
  Group.age Group.gender name age gender height
1        12            F    c  12      F    120
2        14            F    f  14      F    150 #找到了age=14，gender=F这一组中的最大值
3        16            F    e  16      F    160
4        16            M    b  16      M    170
5        18            M    a  18      M    180
by

aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。它无法一次返回若干个统计量。要完成这项任务，可以使用by()函数。格式为：

by(data, INDICES, FUN)

参数 描述
data 是一个数据框或矩阵
INDICES 一个因子或因子组成的列表，定义了分组
FUN 任意函数

例：

构建示例用数据集：

1
2
3
4
5
6
7
8
9
10
11
12
13
> name <- c('a','h','b','c','d','e')
> age <- c(18,16,12,14,16,14)
> gender <- c('M','M',"F",'F',"F",'F')
> height <- c(180,170,120,140,160,150)
> frame1 <- data.frame(name=name,age=age,gender=gender,height=height)
> frame1
  name age gender height
1    a  18      M    180
2    b  16      M    170
3    c  12      F    120
4    d  14      F    140
5    e  16      F    160
6    f  14      F    150

分组，并且使用预设的summar函数进行初步分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> by(frame1, frame1$gender, summary)
frame1$gender: F
     name                age          gender              height     
 Length:4           Min.   :12.0   Length:4           Min.   :120.0  
 Class :character   1st Qu.:13.5   Class :character   1st Qu.:135.0  
 Mode  :character   Median :14.0   Mode  :character   Median :145.0  
                    Mean   :14.0                      Mean   :142.5  
                    3rd Qu.:14.5                      3rd Qu.:152.5  
                    Max.   :16.0                      Max.   :160.0  
--------------------------------------------------------------------------------- 
frame1$gender: M
     name                age          gender              height     
 Length:2           Min.   :16.0   Length:2           Min.   :170.0  
 Class :character   1st Qu.:16.5   Class :character   1st Qu.:172.5  
 Mode  :character   Median :17.0   Mode  :character   Median :175.0  
                    Mean   :17.0                      Mean   :175.0  
                    3rd Qu.:17.5                      3rd Qu.:177.5  
                    Max.   :18.0                      Max.   :180.0

Package：reshape2

reshape2包是一套重构和整合数据集的绝妙的万能工具。

大致说来，你需要首先将数据融合（melt），以使每一行都是唯一的标识符变量组合。然后将数据重铸（cast）为你想要的任何形状。在重铸过程中，你可以使用任何函数对数据进行整合。

例：

构建示例数据：

1
2
3
4
5
6
7
8
9
10
11
12
> library(reshape2)
> ID <- c(1,1,2,3)
> time <- c(1,2,1,2)
> x1 <- c(5,3,6,2)
> x2 <- c(6,5,1,4)
> data1 <- data.frame(ID=ID,time=time,x1=x1,x2=x2)
> data1
  ID time x1 x2
1  1    1  5  6
2  1    2  3  5
3  2    1  6  1
4  3    2  2  4

融合（melt）: 数据集的融合是将它重构为这样一种格式：每个测量变量独占一行，行中带有要唯一确定这个测量所需的标识符变量。

  • 必须指定要唯一确定每个测量所需的变量（ID和time），而表示测量变量名的变量（x1或x2）将由程序为你自动创建。
1
2
3
4
5
6
7
8
9
10
11
> md <- melt(data1, id=c("ID", "time"))
> md
  ID time variable value
1  1    1       x1     5
2  1    2       x1     3
3  2    1       x1     6
4  3    2       x1     2
5  1    1       x2     6
6  1    2       x2     5
7  2    1       x2     1
8  3    2       x2     4

重铸（cast）: dcast()函数读取已融合的数据，并使用你提供的公式和一个（可选的）用于整合数据的函数将其重塑。调用格式为：

newdata <- dcast(md, formula, fun.aggregate)

参数 描述
md 已融合的数据
formula 描述想要的最后结果
fun.aggregate （可选的）数据整合函数。其接受的公式形如：
rowvar1 + rowvar2 + ... ~ colvar1 + colvar2 + ...
在这一公式中，·rowvar1 + rowvar2 + …·定义了要划掉的变量集合，以确定各行的内容，而colvar1 + colvar2 + ...则定义了要划掉的、确定各列内容的变量集合。

这部分的formula在很多R的函数中都会出现，包括回归相关、作图相关等等：

R表达式中常用的符号：

符号 描述
~ 分隔符号，左边为响应变量，右边为解释变量。例如，要通过 x、z 和 w 预测 y，代码为 y ~ x + z + w
+ 分隔预测变量
: 表示预测变量的交互项。例如，要通过 x、z 及 x 与 z 的交互项预测 y，代码为 y ~ x + z + x:z
* 表示所有可能交互项的简洁方式。代码 y~ x * z * w 可展开为 y ~ x + z + w + x:z + x:w + z:w + x:z:w
^ 表示交互项达到某个次数。代码 y ~ (x + z + w)^2 可展开为 y ~ x + z + w + x:z + x:w + z:w
. 表示包含除因变量外的所有变量。例如，若一个数据框包含变量 x、y、z 和 w，代码 y ~ . 可展开为 y ~ x + z + w
- 减号，表示从等式中移除某个变量。例如，y ~ (x + z + w)^2 – x:w 可展开为 y ~ x + z + w + x:z + z:w
-1 删除截距项。例如，表达式 y ~ x - 1 拟合 y 在 x 上的回归，并强制直线通过原点
I() 从算术的角度来解释括号中的元素。例如，y ~ x + (z + w)^2 将展开为 y ~ x + z + w + z:w。
相反, 代码 y ~ x + I((z + w)^2)将展开为 y ~ x + h，h 是一个由 z 和 w 的平方和创建的新变量
function 可以在表达式中用的数学函数。例如，log(y) ~ x + z + w 表示通过 x、z 和 w 来预测 log(y)

例：

image-20200728195014696

2. 频数表和列联表

本节中的数据来自vcd包中的Arthritis数据集。这份数据来自Kock & Edward （1988），表示了一项风湿性关节炎新疗法的双盲临床实验的结果。前几个观测是这样的：

1
2
3
4
5
6
7
8
9
10
> library(vcd)
载入需要的程辑包：grid
> head(Arthritis)
  ID Treatment  Sex Age Improved
1 57   Treated Male  27     Some
2 46   Treated Male  29     None
3 77   Treated Male  30     None
4 17   Treated Male  32   Marked
5 36   Treated Male  46   Marked
6 23   Treated Male  58   Marked
  • 其中的治疗情况（安慰剂治疗、用药治疗）、性别（男性、女性）和改善情况（无改善、一定程度的改善、显著改善）均为类别型因子

生成频数表

R中提供了用于创建频数表和列联表的若干种方法：

函数 描述
table(var1, var2, ..., varN) 使用 N 个类别型变量（因子）创建一个 N 维列联表
xtabs(formula, data) 根据一个公式和一个矩阵或数据框创建一个 N 维列联表
prop.table(table, margins) 依 margins 定义的边际列表将表中条目表示为分数形式
margin is a vector giving the margins to compute sums
for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns.
下同
margin.table(table, margins) 依 margins 定义的边际列表计算表中条目的和
addmargins(table, margins) 将概述边 margins（默认是求和结果）放入表中
ftable(table) 创建一个紧凑的“平铺”式列联表
(1) 一维列联表

可以使用table()函数生成简单的频数统计表。示例如下:

1
2
3
4
5
> data1 <- table(Arthritis$Improved)
> data1

  None   Some Marked 
    42     14     28

可以用prop.table()将这些频数转化为比例值：

1
2
3
4
> prop.table(data1)

     None      Some    Marked 
0.5000000 0.1666667 0.3333333

或使用prop.table()*100转化为百分比：

1
2
3
4
> prop.table(data1)*100

    None     Some   Marked 
50.00000 16.66667 33.33333
(2) 二维列联表
函数table()

对于二维列联表，table()函数的使用格式为：

mytable <- table(A, B)

  • 其中的A是行变量，B是列变量。

例：

构建数据集：

1
2
3
4
5
6
7
8
9
10
11
12
13
> name <- c('a','h','b','c','d','e')
> age <- c(18,16,12,14,16,14)
> gender <- c('M','M',"F",'F',"F",'F')
> height <- c(180,170,120,140,160,150)
> frame1 <- data.frame(name=name,age=age,gender=gender,height=height)
> frame1
  name age gender height
1    a  18      M    180
2    b  16      M    170
3    c  12      F    120
4    d  14      F    140
5    e  16      F    160
6    f  14      F    150

制表：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> table(frame1$name,frame1$age)
   
    12 14 16 18
  a  0  0  0  1
  b  1  0  0  0
  c  0  1  0  0
  d  0  0  1  0
  e  0  1  0  0
  h  0  0  1  0
> table(frame1$name,frame1$gender)
   
    F M
  a 0 1
  b 1 0
  c 1 0
  d 1 0
  e 1 0
  h 0 1

除此之外，xtabs()函数还可使用公式风格的输入创建列联表，格式为：

mytable <- xtabs(~ A + B, data=mydata)

  • 其中的mydata是一个矩阵或数据框。总的来说，要进行交叉分类的变量应出现在公式的右侧（即~符号的右方），以+作为**分隔符**。若某个变量写在公式的左侧，则其为一个频数向量（在数据已经被表格化时很有用）。

例, 对于Arthritis数据，有：

1
2
3
4
5
6
7
8
9
10
> library(vcd)
载入需要的程辑包：grid
> head(Arthritis)
  ID Treatment  Sex Age Improved
1 57   Treated Male  27     Some
2 46   Treated Male  29     None
3 77   Treated Male  30     None
4 17   Treated Male  32   Marked
5 36   Treated Male  46   Marked
6 23   Treated Male  58   Marked
1
2
3
4
5
6
> mytable1 <- xtabs(~ Treatment+Improved, data=Arthritis)
> mytable1
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
1
2
3
4
5
6
> mytable2 <- xtabs(~ Sex+Improved, data=Arthritis)
> mytable2
        Improved
Sex      None Some Marked
  Female   25   12     22
  Male     17    2      6

可以使用margin.table()和prop.table()函数分别生成边际频数和比例。行和与行比例可以这样计算：

  • prop.table(table, margins) & margin.table(table, margins)
  • margin is a vector giving the margins to compute sums for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns.
  • margin除了可以使某个数字之外，也可以直接是某个维度的名称，如：margin.table(mytable1, "Treatment")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
> mytable1
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

# 按照每一行（第1维度）计算比例，即第1维度的比例和为1
> prop.table(mytable1, 1)
         Improved
Treatment      None      Some    Marked
  Placebo 0.6744186 0.1627907 0.1627907
  Treated 0.3170732 0.1707317 0.5121951

# 按照每一列（第2维度）计算比例，即第2维度的比例和为1
> prop.table(mytable1, 2)
         Improved
Treatment      None      Some    Marked
  Placebo 0.6904762 0.5000000 0.2500000
  Treated 0.3095238 0.5000000 0.7500000

# 按照第1，2维度计算，对于多维列联表时适用
> prop.table(mytable1, c(1,2))
         Improved
Treatment None Some Marked
  Placebo    1    1      1
  Treated    1    1      1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> mytable1
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

# 按照每一行（第1维度）求频数和，即每一项都代表原表在第1维度上的频数和
> margin.table(mytable1, 1)
Treatment
Placebo Treated 
     43      41 

# 按照每一列（第2维度）求频数和，即每一项都代表原表在第2维度上的频数和
> margin.table(mytable1, 2)
Improved
  None   Some Marked 
    42     14     28 

# 按照第1，2维度计算，对于多维列联表时适用
> margin.table(mytable1, c(1,2))
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

可以使用addmargins()函数为表格添加边际和（即求每一行每一列的和，并表示在原表的边界上）。

例如，以下代码添加了各行的和与各列的和：

1
2
3
4
5
6
7
8
9
10
11
12
> mytable1
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

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

在使用addmargins()时，默认行为是为表中所有的变量创建边际和。可以通过设置参数进行对所有行或者所有列进行求和：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 对每一列进行求和
> addmargins(mytable1, 1)
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
  Sum       42   14     28

# 对每一行进行求和
> addmargins(mytable1, 2)
         Improved
Treatment None Some Marked Sum
  Placebo   29    7      7  43
  Treated   13    7     21  41
Package：gmodels

使用gmodels包中的CrossTable()函数也是创建二维列联表的方法。

  • CrossTable()函数仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二维列联表。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
> library(gmodels)
Warning message:
程辑包‘gmodels’是用R版本4.0.2 来建造的 
> CrossTable(Arthritis$Treatment, Arthritis$Improved)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  84 

 
                    | Arthritis$Improved 
Arthritis$Treatment |      None |      Some |    Marked | Row Total | 
--------------------|-----------|-----------|-----------|-----------|
            Placebo |        29 |         7 |         7 |        43 | 
                    |     2.616 |     0.004 |     3.752 |           | 
                    |     0.674 |     0.163 |     0.163 |     0.512 | 
                    |     0.690 |     0.500 |     0.250 |           | 
                    |     0.345 |     0.083 |     0.083 |           | 
--------------------|-----------|-----------|-----------|-----------|
            Treated |        13 |         7 |        21 |        41 | 
                    |     2.744 |     0.004 |     3.935 |           | 
                    |     0.317 |     0.171 |     0.512 |     0.488 | 
                    |     0.310 |     0.500 |     0.750 |           | 
                    |     0.155 |     0.083 |     0.250 |           | 
--------------------|-----------|-----------|-----------|-----------|
       Column Total |        42 |        14 |        28 |        84 | 
                    |     0.500 |     0.167 |     0.333 |           | 
--------------------|-----------|-----------|-----------|-----------|

CrossTable()函数有很多选项，可以做许多事情：计算（行、列、单元格）的百分比；指定小数位数；进行卡方、Fisher和McNemar独立性检验；计算期望和（皮尔逊、标准化、调整的标准化）残差；将缺失值作为一种有效值；进行行和列标题的标注；生成SAS或SPSS风格的输出。

这里不赘述，因为不常用，具体可参阅help(CrossTable)

(3) 多维列联表

如果有两个以上的类别型变量，那么就是在处理多维列联表。之前提到的table()和xtabs()函数同样可以用来生成多维列联表，只需要改变公式即可：

使用table()创建时，不同的维度需要通过,隔开：table(var1, var2, var3,...)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
> attach(Arthritis)

# 这里的输入变量的顺序决定了输出表格的不同维度：
# 第一维度：Treatment；第二维度：Sex，...
> table(Treatment,Sex,Improved)
, , Improved = None

         Sex
Treatment Female Male
  Placebo     19   10
  Treated      6    7

, , Improved = Some

         Sex
Treatment Female Male
  Placebo      7    0
  Treated      5    2

, , Improved = Marked

         Sex
Treatment Female Male
  Placebo      6    1
  Treated     16    5

> detach(Arthritis)

通过xtab()创建时，不同的变量需要通过+隔开：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 这里的输入变量的顺序决定了输出表格的不同维度：
# 第一维度：Treatment；第二维度：Sex，...
> xtabs(~ Treatment+Sex+Improved, data=Arthritis)
, , Improved = None

         Sex
Treatment Female Male
  Placebo     19   10
  Treated      6    7

, , Improved = Some

         Sex
Treatment Female Male
  Placebo      7    0
  Treated      5    2

, , Improved = Marked

         Sex
Treatment Female Male
  Placebo      6    1
  Treated     16    5

此时，ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表。

1
2
3
4
5
6
7
8
> table1 <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
> ftable(table1)
                 Improved None Some Marked
Treatment Sex                             
Placebo   Female            19    7      6
          Male              10    0      1
Treated   Female             6    5     16
          Male               7    2      5

在多维列联表的状态下，边界求和函数addmargins()，margin.table()，prop.table()等仍然适用：

  • 在多维列联表的情况下，margin的使用更加自由:
    • 当margin为某个数字的时候，指的是对某个维度进行求和
    • 当margin为某个向量的时候，可以是通过数字指定某个维度，如``

margin.table()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
> table1 <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
# 对第一维度（本例中为Treatment）求和
> margin.table(table1, 1)
Treatment
Placebo Treated 
     43      41 

> margin.table(table1, "Treatment")
Treatment
Placebo Treated 
     43      41

# 对第二维度（本例中为Sex）求和
> margin.table(table1, 2)
Sex
Female   Male 
    59     25 

# 对第三维度（本例中为Improved）求和
> margin.table(table1, 3)
Improved
  None   Some Marked 
    42     14     28 

# 对1，2维度求和，即重新构建次级表格
> margin.table(table1, c(1,2))
         Sex
Treatment Female Male
  Placebo     32   11
  Treated     27   14

# 对2，3维度求和，不同的顺序可以改变输出表格的行列
> margin.table(table1, c(2,3))
        Improved
Sex      None Some Marked
  Female   25   12     22
  Male     17    2      6

> margin.table(table1, c(3,2))
        Sex
Improved Female Male
  None       25   17
  Some       12    2
  Marked     22    6

prob.table(): 为了方便展示，这里使用了ftable()

  • margin使用的方法同margin.table
1
2
3
4
5
6
7
8
# 按照第一维度进行求概率
> ftable(prop.table(table1, 1))
                 Improved       None       Some     Marked
Treatment Sex                                             
Placebo   Female          0.44186047 0.16279070 0.13953488
          Male            0.23255814 0.00000000 0.02325581
Treated   Female          0.14634146 0.12195122 0.39024390
          Male            0.17073171 0.04878049 0.12195122

addmargins():

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 计算原表1，2维度上各项的概率，同时计算各个项在第3维度的概率和
> ftable(addmargins(prop.table(table1, c(1, 2)), 3))
                 Improved       None       Some     Marked        Sum
Treatment Sex                                                        
Placebo   Female          0.59375000 0.21875000 0.18750000 1.00000000
          Male            0.90909091 0.00000000 0.09090909 1.00000000
Treated   Female          0.22222222 0.18518519 0.59259259 1.00000000
          Male            0.50000000 0.14285714 0.35714286 1.00000000

# 计算原表1，2维度上各项的概率，同时计算各个项在第2维度的概率和
> ftable(addmargins(prop.table(table1, c(1, 2)), 2))
                 Improved       None       Some     Marked
Treatment Sex                                             
Placebo   Female          0.59375000 0.21875000 0.18750000
          Male            0.90909091 0.00000000 0.09090909
          Sum             1.50284091 0.21875000 0.27840909
Treated   Female          0.22222222 0.18518519 0.59259259
          Male            0.50000000 0.14285714 0.35714286
          Sum             0.72222222 0.32804233 0.94973545

二、单项描述性统计分析

数学函数

函数 描述
abs(x) 绝对值
sqrt() 平方根
ceiling(x) 不小于 x 的最小整数, ceiling(3.475)返回值为 4
floor(x) 不大于 x 的最大整数, floor(3.475)返回值为 3
trunc(x) 向 0 的方向截取的 x 中的整数部分, trunc(5.99)返回值为 5
round(x, digits=n) 将 x 舍入为指定位的小数, round(3.475, digits=2)返回值为 3.48
signif(x, digits=n) 将 x 舍入为指定的有效数字位数, signif(3.475, digits=2)返回值为 3.5
cos(x)、sin(x)、tan(x) 余弦、正弦和正切
acos(x)、asin(x)、atan(x) 反余弦、反正弦和反正切
cosh(x)、sinh(x)、tanh(x) 双曲余弦、双曲正弦和双曲正切
acosh(x)、asinh(x)、atanh(x) 反双曲余弦、反双曲正弦和反双曲正切
log(x,base=n)
log(x)
log10(x)		 对 x 取以 n 为底的对数
为了方便起见：
log(x)为自然对数
log10(x)为常用对数
log(10)返回值为 2.3026
log10(10)返回值为 1
exp(x) 指数函数

统计函数

函数 描述
mean(x) 平均数，x为数值向量
median(x) 中位数
sd(x) 样本标准差: $\sum / (n-1)$
这里的Standard Deviation是对样本的标准差，即Sample Standard Deviation，下面说的Variance也是一样
var(x) 样本方差: 由Central Limit Theorem, $\sigma/\sqrt{n}$
$\sigma$是总体方差
mad(x) 绝对中位差（median absolute deviation）
quantile(x,probs) 求分位数。其中 x 为待求分位数的数值型向量，probs 为一个由[0,1]之间的概率值组成 的数值向量
求 x 的 30%和 84%分位点:
y <- quantile(x, c(.3,.84))
range(x) 求值域
sum(x) 求和
diff(x, lag=n) 滞后差分，lag 用以指定滞后几项。默认的 lag 值为 1
x <- c(1, 5, 23, 29)
diff(x)返回值为 c(4, 18, 6)
min(x) 求最小值
max(x) 求最大值
scale(x,center=TRUE, scale=TRUE) 为数据对象 x 按列进行中心化(center=TRUE)或标准化(center=TRUE,scale=TRUE)
软件使用总结 > R
#速查 #R #Statistics
03. R <- 基本统计分析Ⅰ - 简要分析、频数表和列联表、基本数学和统计函数
https://zhenyumi.github.io/posts/c5486dc5/
作者
向海
发布于
2020年7月25日
更新于
2020年8月1日
许可协议