描述性统计
梁雪枫
统计分析分为统计描述和统计推断两个部分，统计描述是通过绘制统计图、计算统计量等方法描述数据的分布特征，是数据分析的基本步骤。
常用统计量
矩
设X和Y是随机变量，若\[E(X^{k}),k=1,2,\cdot \cdot \cdot \]存在，则称它为X的k阶原点矩，简称k阶矩。 若\[E\left \{ [X-E(x)]^{k} \right \},k=2,3,\cdot \cdot \cdot\]存在，则称它为X的k阶中心距。 若\[E(X^{k}Y^{l}),k,l=1,2,\cdot \cdot \cdot\]存在，则称它为X和Y的k+l阶混合距。 若\[E\left \{ [X-E(X)]^{k}[Y-E(Y)]^{l} \right \},k,l=1,2,\cdot \cdot \cdot\]存在，则称它为X和Y的k+l阶混合中心距。 X的数学期望E(X)是X的一阶原点矩，方差D(X)是X的二阶中心矩，协方差Cov(X,Y)是X和Y的二阶混合中心矩。 ### 均值(Mean) 一阶原点矩又称均数是一组数据的平均值,均数（记为\(\bar{x}\)）定义为 \[\bar{X}=\frac{1}{n}\sum _{i=1}^{n}X_{i}\] 用它来描述正态分布数据的集中趋势。 ### 标准差(Standard Deviation) 样本方差定义为\[S^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}\],标准差是方差的算术平方根，是一组数值自均数分散开来的程度的一种测量观念。定义为 \[\delta = \sqrt{\frac{1}{N}\sum_{i=1}^{n}(x_i-u)^2}\] 一个较大的标准差，代表大部分的数值和其平均值之间差异较大；一个较小的标准差，代表这些数值较接近平均值。 例 已知50名患者的收缩压（mmHg）分别为： 147 163 159 124 120 94 135 185 109 143 116 129 157 146 149 127 124 160 101 129 130 154 151 119 128 147 127 122 145 159 141 131 117 139 142 152 147 157 134 146 144 119 160 136 122 172 170 109 151 144 求血压的集中趋势和离散趋势及集中趋势的95%可信区间。 思路：先判断数据是否为正态分布，然后根据结果选择描述集中趋势的统计量。
sbp <- c(147,163,159,124,120,94,135,185,109,143,116,129,157,146,149,127,124,160,101,129,130,154,151,119,128,147,127,122,145,159,141,131,117,139,142,152,147,157,134,146,144,119,160,136 ,122,172,170,109,151,144)
qqnorm(sbp)
qqline(sbp)
#plot(density(sbp)) 核密度
hist(sbp,freq = F) #freq=T,则绘制频数
result <- shapiro.test(sbp)
result##
## Shapiro-Wilk normality test
##
## data: sbp
## W = 0.99175, p-value = 0.9783ks.test(sbp,"pnorm",mean(sbp),sd(sbp))## Warning in ks.test(sbp, "pnorm", mean(sbp), sd(sbp)): ties should not be
## present for the Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: sbp
## D = 0.071566, p-value = 0.9599
## alternative hypothesis: two-sided样本大小在3和5000之间,选择Shapiro-Wilk进行正态性检验。W值为0.9917495，P值为0.9782552大于0.05，不能拒绝其于正态分布一致的假设。如果样本数较大，可以选择Kolmogorov-Smirnov检验。对于服从正态分布的数据，选择均数和标准差描述其集中趋势和离散趋势。
mean(sbp)## [1] 138.64sd(sbp)## [1] 18.91308均值为138.64,标准差为18.9130772。对均数和标准差的计算还可以通过base包中的summary()函数，该函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。epicalc包中的summ()和mosaic包中的favstats()等函数也可获得类似结果，如favstats()一次就可以完成均数和标准差的计算。
favstats(sbp)## min Q1 median Q3 max mean sd n missing
## 94 124.75 141.5 151 185 138.64 18.91308 50 0对均数的95%可信区间的计算可通过t.test()获得，对一个给定的可信区间，它表示一个总体参数的估计范围。根据中位数和均数可以快速检查数据分布，如中位数小于均数，说明分布有可能向右倾斜。
t.test(sbp)##
## One Sample t-test
##
## data: sbp
## t = 51.834, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 133.265 144.015
## sample estimates:
## mean of x
## 138.64通过设置conf.level=0.99，可以将可信区间水平提高到99%。
t.test(sbp,conf.level=0.99)##
## One Sample t-test
##
## data: sbp
## t = 51.834, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
## 131.4719 145.8081
## sample estimates:
## mean of x
## 138.64中位数(Median)
对于一组有限个数的数据来说，它们的中位数是这样的一种数：这群数据里的一半的数据比它大，而另外一半数据比它小。 计算有限个数的数据的中位数的方法是：把所有的同类数据按照大小的顺序排列。如果数据的个数是奇数，则中间那个数据就是这群数据的中位数；如果数据的个数是偶数，则中间那2个数据的算术平均值就是这群数据的中位数。通常用中位数来描述非正态分布数据的集中趋势，极值对中位数影响不大。定义为 实数\(x_1\), \(x_2\), \(\dots\), \(x_n\)按大小顺序（顺序，降序皆可）排列为\(x'_1\), \(x'_2\), \(\dots\) , \(x'_n\)、实数数列\(x=(x_1, x_2, \dots , x_n)\)的中位数\(\mathrm{Q}_\frac{1}{2}(x)\)为 \[\mathrm{Q}_\frac{1}{2}(x) = \begin{cases} x'_\frac{n + 1}{2}, & \mbox{if } n \mbox{ is odd.} \\ \frac{1}{2}( x'_\frac{n}{2} + x'_{\frac{n}{2} + 1}), & \mbox{if } n \mbox{ is even.} \end{cases} \] odd 为奇数，even 为偶数。 ###四分位差(quartile deviation) 是上四分位数（QU，即位于75%）与下四分位数（QL，即位于25%）的差的一半。计算公式为：\(Qd=QU-QL\)。四分位差反映了中间50%数据的离散程度，其数值越小，说明中间的数据越集中；其数值越大，说明中间的数据越分散。四分位差不受极值的影响。 例 某地17名患者的月收入分别为： 23408 3468 1939 4360 23545 12233 4583 3546 35781 6578 8981 1345 5567 23455 23564 7623 14334 求收入的集中趋势和离散趋势及集中趋势的95%可信区间。 思路：仍然先判断数据是否为正态分布，然后根据结果选择描述集中趋势的统计量。
income <- c(23408,3468,1939,4360,23545,12233,4583,3546,35781,6578,8981,1345,5567,
23455, 23564,7623,14334)
qqnorm(income)
qqline(income)
#plot(density(income))
hist(income,freq = F)
result <- shapiro.test(income)
result##
## Shapiro-Wilk normality test
##
## data: income
## W = 0.85488, p-value = 0.01274根据QQ图和Shapiro-Wilk进行正态性检验结果，W值为0.8548814，P值为0.0127409小于0.05，选择中位数和四分位差描述其集中趋势和离散趋势。
quantile(income)## 0% 25% 50% 75% 100%
## 1345 4360 7623 23408 35781中位数为7623，四分位差为75%的分位数减去25%的分位数。对中位数的可信区间估计，可通过wilcox.test()函数获得
wilcox.test(income, conf.int=TRUE)##
## Wilcoxon signed rank test
##
## data: income
## V = 153, p-value = 1.526e-05
## alternative hypothesis: true location is not equal to 0
## 95 percent confidence interval:
## 5460.0 17820.5
## sample estimates:
## (pseudo)median
## 12376.5数学期望(mathematical expectation)
离散型随机变量：离散型随机变量的一切可能的取值\(x_{i}\)与对应的概率\(P_{i}(=x_{i})\)之积的和称为该离散型随机变量的数学期望，记为\(E(X)\)。数学期望是最基本的数学特征之一。它反映随机变量平均取值的大小 \[E(X)=\sum_{i}x_{i}p_{i}\] 连续型随机变量：若随机变量X的分布函数\(F(x)\)可表示成一个非负可积函数\(f(x)\)的积分，则称\(X\)为连续性随机变量，\(f(x)\)称为\(X\)的概率密度函数，积分值为X的数学期望，记为\(E(X)\)。 \[E(X)=\int_{-\infty }^{+\infty }xf(x)dx\] ###方差(Variance) 方差是各个数据与平均数之差的平方的平均数。在概率论和数理统计中，方差用来度量随机变量和其数学期望（即均值）之间的偏离程度。 设\(X\)为随机变量，如果\(E{[X-E(X)]^2}\)存在，则称\(E{[X-E(X)]^2}\)为\(X\)的方差，记为\(Var(X)\)。 离散型随机变量方差计算公式为\[Var(X)=E(X^{2})-(E(X))^{2}\],连续型随机变量方差计算公式为\(Var(x)=\int_{-\infty}^{+\infty} (x-E(X))^{2}f(x) \text{d}x=E(X^{2})-(E(X))^{2}\) 例 计算样本(2,5,78,45,89,124)的方差
s <- c(2,5,78,45,89,124)
var(s)## [1] 2365.367众数(Mode)
观察资料中出现次数最多的数值或类别，不受极值影响，由于可能有不只一个，也可能没有众数，一般不适合进行统计分析。 例 计算样本(4,22,31,33,3,27,27,27,27,569,110,8,21,31,33,33)的众数
S <- c(4,22,31,33,3,27,27,27,27,569,110,8,21,31,33,33)
names(which.max(table(S)))## [1] "27"协方差(Covariance)
协方差用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。设\(X\),\(Y\)为两个随机变量，称\(E{[X-E(X)][Y-E(Y)]}\)为\(X\)和\(Y\)的协方差，记录\(Cov(X,Y)\)。方差是协方差的一种特殊情况，即当两个变量是相同的情况。 \[Cov(X,Y)=E\left \{ [X-E(X)][Y-E(Y)] \right \}=E(XY)-E(X)E(Y)\] 例 计算X(2,5,7)和Y(6,7,9)的协方差。
x <- c(2,5,7)
y <- c(6,7,9)
cov(x,y)## [1] 3.666667相关系数(Correlation coefficient)
相关系数是用以反映变量之间相关关系密切程度的统计指标。相关系数是按积差方法计算，同样以两变量与各自平均值的离差为基础，通过两个离差相乘来反映两变量之间相关程度。当\(Var(X)>0, Var(Y)>0\)时，称\(Cov(X,Y)/sqrt(Var(X)*Var(Y))\)为\(X\)与\(Y\)的相关系统。 \[\rho (X,Y)=\frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}}\] 例 计算X(2,5,7)和Y(6,7,9)的相关系数
x <- c(2,5,7)
y <- c(6,7,9)
cor(x,y)## [1] 0.953821偏度(skewness)
是统计数据分布偏斜方向和程度的度量，是统计数据分布非对称程度的数字特征。设分布函数\(F(x)\)有中心矩\(\mu_{2}=E(X-E(X))^2\),\(\mu_{3}=E(X-E(X))^3\),则\(\frac{\mu_{3}}{\mu_{2}^{\frac{3}{2}}}\)为偏度系数。当\(C_{s}\)>0时，概率分布偏向均值右则,\(C_{s}\)<0时，概率分布偏向均值左则。 \[C_{s}=\frac{\mu_{3}}{\mu_{2}^{\frac{3}{2}}}\] ###峰度(kurtosis) 表征概率密度分布曲线在平均值处峰值高低的特征数。峰度刻划不同类型的分布的集中和分散程序。设分布函数\(F(x)\)有中心矩\(\mu_{2}=E(X-E(X))^2\), \(\mu_{4}=E(X-E(X))^4\)，则\(C_{k}=\frac{\mu_{4}}{\mu_{2}^{2}}-3\)为峰度系数。 \[C_{k}=\frac{\mu_{4}}{\mu_{2}^{2}}-3\] 例 计算10000个正态分布的样本的偏度和峰度
S<-rnorm(10000)
skewness(S)## [1] -0.01502003kurtosis(S)## [1] -0.06567366hist(S,breaks=100)
###几何平均数（Geometric mean） 是n个变量值连乘积的n次方根,是用于反映一组经对数转换后呈对称分布的变量值在数量上的平均水平即对数正态分布数据，在医学研究中常适用于免疫学的指标。对于变量值呈倍数关系或呈对数正态分布（正偏态分布），如抗体效价及抗体滴度，某些传染病的潜伏期，细菌计数等，宜用几何均数表示其平均水平。 \[H=G=\sqrt[n]{X_{1}*X_{2}*...*X_{n}}=\sum \sqrt[n]{\prod_{i=1}^{n}X_{n}}\] 例 5名学龄儿童的麻疹血凝抑制抗体滴度为1：25，1：50，1：50，1：100，1：400，求几何均数及标准差。
geomean <- function(x, na.rm = FALSE, trim = 0, ...)
{
exp(mean(log(x, ...), na.rm = na.rm, trim = trim, ...))
}
geosd <- function(x, na.rm = FALSE, ...)
{
exp(sd(log(x, ...), na.rm = na.rm, ...))
}
s<-c(25,50,50,100,400)
shapiro.test(log(s))##
## Shapiro-Wilk normality test
##
## data: log(s)
## W = 0.91408, p-value = 0.4925geomean(s)## [1] 75.78583geosd(s)## [1] 2.86111也可以安装NCStats包，调用geomean和geosd()函数。 ###变异系数（Coefficient of Variation） 是刻画数据相对分散性的一种度量，记为\(c_v\),是概率分布离散程度的一个归一化量度，其定义为标准差\(\ \sigma\)与平均值\(\ \mu\) 之比 \[c_v = {\sigma \over \mu }\] ###样本校正平方和（CSS） 样本与均值差的平方的求和\(CSS=\sum_{i=1}^{n}(x_{i}-\overline{x})\) ###样本未校正平方和（USS） 样本值平方的求和\(USS=\sum_{i=1}^{n}x_{i}^{2}\) ###标准误（Standard Deviation） 是某种统计量在抽样分布上的标准差称为该种统计量的标准误，即样本统计量的标准差，是描述对应的样本统计量抽样分布的离散程度及衡量对应样本统计量抽样误差大小的尺度。设n个测量值的误差为v1、v2……vn，则这组测量值的标准误差\(\sigma\) \[\sigma =\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}v_{i}^2}\] ###极差 描述样本分散性，数据越分散，其极差越大。 \[R=max(x)-min(x)\] 例 对例1求变异系数、样本校正平方和、样本未校正平方和、极差和均值的标准误。
cv <- 100*sd(sbp)/mean(sbp)
cv## [1] 13.64186css <- sum((sbp-mean(sbp))^2)
css## [1] 17527.52uss <- sum(sbp^2)
uss## [1] 978580r <- max(sbp)-min(sbp)
r## [1] 91通常由于总体的均数或总体的方差并不知道，样本均值的标准误\(SD=\frac{s}{\sqrt{n}}\),s为标准差，n为样本数。
sd <- sd(sbp)/sqrt(length(sbp))
sd## [1] 2.674713数据中心化和标准化
数据中心化是将某变量中的观察值减去该变量的平均数，数据标准化将某变量中的观察值减去该变量的平均数，然后除以该变量的标准差。经标准化的数据都是没有单位的纯数量。对变量进行的标准差标准化可以消除量纲（单位）影响和变量自身变异的影响。 例 对下表中三科成绩进行标准化。 Math Science English — — — 502 95 25 465 67 12 621 78 22 575 66 18 454 96 15 634 89 30 576 78 37 421 56 12 599 68 22 666 100 38
R语言中scale()函数可以实现数据标准化，两个参数center和scale为True分别表示计算中心化和标准化
Math <- c(502,465,621,575,454,634,576,421,599,666)
Science <- c(95,67,78,66,96,89,78,56,68,100)
English <- c(25,12,22,18,15,30,37,12,22,38)
Student <- as.data.frame(cbind(Math,Science,English))
options(digits=2) #限定为2位小数
scale(Student[,1:3],center = T,scale = F) #数据中心化 ## Math Science English
## [1,] -49 15.7 1.9
## [2,] -86 -12.3 -11.1
## [3,] 70 -1.3 -1.1
## [4,] 24 -13.3 -5.1
## [5,] -97 16.7 -8.1
## [6,] 83 9.7 6.9
## [7,] 25 -1.3 13.9
## [8,] -130 -23.3 -11.1
## [9,] 48 -11.3 -1.1
## [10,] 115 20.7 14.9
## attr(,"scaled:center")
## Math Science English
## 551 79 23scale(Student[,1:3],center = F,scale = T) #数据标准化## Math Science English
## [1,] 0.85 1.12 0.96
## [2,] 0.79 0.79 0.46
## [3,] 1.06 0.92 0.84
## [4,] 0.98 0.78 0.69
## [5,] 0.77 1.13 0.57
## [6,] 1.08 1.05 1.15
## [7,] 0.98 0.92 1.42
## [8,] 0.72 0.66 0.46
## [9,] 1.02 0.80 0.84
## [10,] 1.13 1.18 1.45
## attr(,"scaled:scale")
## Math Science English
## 587 85 26apply()函数或sapply()函数计算所选择的任意描述性统计量。对于sapply()函数，其使用格式为：sapply(x,FUN,options)其中的x是输入的数据框（或矩阵），FUN为一个任意的函数。如果指定了options，它们将被传递给FUN。你可以在这里插入的典型函数有mean、sd、var、min、max、median、length、range和quantile。可以根据需要自定义需要的统计量，如下
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))
}
data(drugDat,package = "elrm")
sapply(drugDat,mystats)## sex treatment recovered n
## n 4.00 4.00 4.0 4.00
## mean 0.50 0.50 11.5 24.75
## stdev 0.58 0.58 3.9 5.91
## skew 0.00 0.00 0.0 0.18
## kurtosis -2.44 -2.44 -2.1 -2.16Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。pastecs包中有一个名为stat.desc()的函数，它可以计算种类繁多的描述性统计量。使用格式为：stat.desc(x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95)其中的x是一个数据框或时间序列。若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。若desc=TRUE（同样也是默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。最后，若norm=TRUE（不是默认的），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro–Wilk正态检验结果。这里使用了p值来计算平均数的置信区间（默认置信度为0.95）。psych包也拥有一个名为describe()的函数，它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误。
数据操作
数据输入
readr包中的函数使数据读入的速度更快，相对于基础包中的函数，对字符类型并不需要指定stringsAsFactors = FALSE防止字符类型自动转为因子，对列名限制更少。固定分割的数据使用read_delim(), read_csv(), read_tsv()和read_csv2()函数，固定宽度的数据使用read_fwf()和 read_table()。
WHO<- read_csv("WHO.csv",col_names=T) #col_names相当于header=T，默认为True
#可以从压缩包或网站上直接输入
mtcars <- read_csv(system.file("extdata/mtcars.csv.bz2", package = "readr"))
#mtcars <- read_csv("https://github.com/hadley/readr/raw/master/inst/extdata/mtcars.csv")数据输出
用readr包读入的数据，变量的引用使用如下格式WHO$Adolescent fertility rate (%)，不同于通常的引用。write_csv()将数据框快速的输出为csv文件。
who <- read_csv("WHO.csv", col_types = list(
CountryID = col_integer(),
Continent=col_double(),
Country=col_factor(c("Country")) #col_date() 使用Y-m-d格式，col_datetime()使用 ISO8601日期时间格式
))## Warning: 202 parsing failures.
## row col expected actual
## 1 Country value in level set Afghanistan
## 2 Country value in level set Albania
## 3 Country value in level set Algeria
## 4 Country value in level set Andorra
## 5 Country value in level set Angola
## ... ....... .................. ...........
## .See problems(...) for more details.class(who) #tbl_df、tbl和data.frame类型## [1] "tbl_df" "tbl" "data.frame"字符串操作
合并字符串
IT <- c("google","baidu","bing")
res <- str_c(1:3,IT,sep=' ',collapse=' ')
str_c('My work place is ',res,collapse=' ')## [1] "My work place is 1 google 2 baidu 3 bing"计算字符串长度
str_length(c("programming R and Python", 123,res))## [1] 24 3 23按位置取子字符串
str_sub(IT, 1, 3)## [1] "goo" "bai" "bin"子字符串重新赋值
capital <-toupper(str_sub(IT,1,1))
str_sub(IT, rep(1,3),rep(1,3)) <- capital 重复字符串
str_dup(IT, c(2,3,4))## [1] "GoogleGoogle" "BaiduBaiduBaidu" "BingBingBingBing"加空白和去除空白
str_pad(IT, 10, "both")## [1] " Google " " Baidu " " Bing "str_trim(IT)## [1] "Google" "Baidu" "Bing"根据正则表达式检验是否匹配
str_detect(IT, "g$") #查找以g结尾## [1] FALSE FALSE TRUEstr_detect(IT, "[aiu]") #查找是否包含a、i、u## [1] FALSE TRUE TRUE查找匹配的字符串位置
str_locate(IT, "a") #返回起始和结束的位置## start end
## [1,] NA NA
## [2,] 2 2
## [3,] NA NA提取匹配的部分
str_extract(IT, "[a-z]+")## [1] "oogle" "aidu" "ing"str_extract(IT, "[a-z]{1,3}")## [1] "oog" "aid" "ing"str_match(IT, "[a-z]+")## [,1]
## [1,] "oogle"
## [2,] "aidu"
## [3,] "ing"替换匹配的部分
str_replace(IT, "[aeiou]", "-")## [1] "G-ogle" "B-idu" "B-ng"分割
str_split(res, " ")## [[1]]
## [1] "1" "google" "2" "baidu" "3" "bing"数据操作
dplyr包将plyr包中的ddply()等函数进一步分离强化,专注接受dataframe对象,大幅提高了运算速度,并且提供了更稳健的与其它数据库对象间的接口。 ####数据集类型 将过长过大的数据集转换为显示更友好的tbl_df类型
iris_df<- tbl_df(iris)
iris_df## Source: local data frame [150 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 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
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## .. ... ... ... ... ...筛选
filter 用于选择满足条件的观测（行），第一个参数是 data frame 名字，第二个参数是条件。
filter(iris_df, Species == "versicolor") #选取 Species == versicolor的观测## Source: local data frame [50 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 1 7.0 3.2 4.7 1.4 versicolor
## 2 6.4 3.2 4.5 1.5 versicolor
## 3 6.9 3.1 4.9 1.5 versicolor
## 4 5.5 2.3 4.0 1.3 versicolor
## 5 6.5 2.8 4.6 1.5 versicolor
## 6 5.7 2.8 4.5 1.3 versicolor
## 7 6.3 3.3 4.7 1.6 versicolor
## 8 4.9 2.4 3.3 1.0 versicolor
## 9 6.6 2.9 4.6 1.3 versicolor
## 10 5.2 2.7 3.9 1.4 versicolor
## .. ... ... ... ... ...filter(iris_df, Sepal.Length %in% c(7.0, 5.2,6.6)) #选取Sepal.Length为7.0，5.2，6.6的观测## Source: local data frame [7 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 1 5.2 3.5 1.5 0.2 setosa
## 2 5.2 3.4 1.4 0.2 setosa
## 3 5.2 4.1 1.5 0.1 setosa
## 4 7.0 3.2 4.7 1.4 versicolor
## 5 6.6 2.9 4.6 1.3 versicolor
## 6 5.2 2.7 3.9 1.4 versicolor
## 7 6.6 3.0 4.4 1.4 versicolor对于多条件的选择，需要完整条件的，然后使用集合运算符将条件拼接起来。集合运算符有 !、|、&、xor(交补)。条件的判断符有>(=)、<(=)、==、!=、%in% (判断元素是否在集合或者列表内，返回逻辑值)。
filter(iris_df, Sepal.Length>=6.3 & Species=="versicolor")## Source: local data frame [14 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 1 7.0 3.2 4.7 1.4 versicolor
## 2 6.4 3.2 4.5 1.5 versicolor
## 3 6.9 3.1 4.9 1.5 versicolor
## 4 6.5 2.8 4.6 1.5 versicolor
## 5 6.3 3.3 4.7 1.6 versicolor
## 6 6.6 2.9 4.6 1.3 versicolor
## 7 6.7 3.1 4.4 1.4 versicolor
## 8 6.3 2.5 4.9 1.5 versicolor
## 9 6.4 2.9 4.3 1.3 versicolor
## 10 6.6 3.0 4.4 1.4 versicolor
## 11 6.8 2.8 4.8 1.4 versicolor
## 12 6.7 3.0 5.0 1.7 versicolor
## 13 6.7 3.1 4.7 1.5 versicolor
## 14 6.3 2.3 4.4 1.3 versicolor排列
arrange 用于根据变量排序，如果排序依据（列）是字符，按照字母表的顺序，如果是数字，默认按照从小到大的顺序排序，如果需要使用逆序排，可以使用desc(var) 或者 -var。
arrange(iris_df, Petal.Length)## Source: local data frame [150 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 1 4.6 3.6 1.0 0.2 setosa
## 2 4.3 3.0 1.1 0.1 setosa
## 3 5.8 4.0 1.2 0.2 setosa
## 4 5.0 3.2 1.2 0.2 setosa
## 5 4.7 3.2 1.3 0.2 setosa
## 6 5.4 3.9 1.3 0.4 setosa
## 7 5.5 3.5 1.3 0.2 setosa
## 8 4.4 3.0 1.3 0.2 setosa
## 9 5.0 3.5 1.3 0.3 setosa
## 10 4.5 2.3 1.3 0.3 setosa
## .. ... ... ... ... ...arrange(iris_df, desc(Petal.Length))## Source: local data frame [150 x 5]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## (dbl) (dbl) (dbl) (dbl) (fctr)
## 1 7.7 2.6 6.9 2.3 virginica
## 2 7.7 3.8 6.7 2.2 virginica
## 3 7.7 2.8 6.7 2.0 virginica
## 4 7.6 3.0 6.6 2.1 virginica
## 5 7.9 3.8 6.4 2.0 virginica
## 6 7.3 2.9 6.3 1.8 virginica
## 7 7.2 3.6 6.1 2.5 virginica
## 8 7.4 2.8 6.1 1.9 virginica
## 9 7.7 3.0 6.1 2.3 virginica
## 10 6.3 3.3 6.0 2.5 virginica
## .. ... ... ... ... ...选择
select 用于选择列,类似于R自带的 subset() 函数,select中负号表示不选择。其中变量的声明还有其他形式，比如B:F表示从 B 列到 F 列所有列；ends_with(“string”) 表示选取列名以 string 结尾的全部列；contains(“string”) 表示选取列名中含有 string 的所有列。
select(iris_df, Petal.Length)## Source: local data frame [150 x 1]
##
## Petal.Length
## (dbl)
## 1 1.4
## 2 1.4
## 3 1.3
## 4 1.5
## 5 1.4
## 6 1.7
## 7 1.4
## 8 1.5
## 9 1.4
## 10 1.5
## .. ...变形
mutate用于添加新的变量，直接使用列名进行计算得到新变量即可。可使用刚添加的变量，也就是在一个语句中可以多个变量，而且变量可以来源于刚新建的变量。
mutate(iris_df, double=Petal.Length*2,quadruple=double*2)## Source: local data frame [150 x 7]
##
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species double
## (dbl) (dbl) (dbl) (dbl) (fctr) (dbl)
## 1 5.1 3.5 1.4 0.2 setosa 2.8
## 2 4.9 3.0 1.4 0.2 setosa 2.8
## 3 4.7 3.2 1.3 0.2 setosa 2.6
## 4 4.6 3.1 1.5 0.2 setosa 3.0
## 5 5.0 3.6 1.4 0.2 setosa 2.8
## 6 5.4 3.9 1.7 0.4 setosa 3.4
## 7 4.6 3.4 1.4 0.3 setosa 2.8
## 8 5.0 3.4 1.5 0.2 setosa 3.0
## 9 4.4 2.9 1.4 0.2 setosa 2.8
## 10 4.9 3.1 1.5 0.1 setosa 3.0
## .. ... ... ... ... ... ...
## Variables not shown: quadruple (dbl)分类汇总
summarise可以用于分类汇总,实际上它是把 data frame 依据分组依据拆分成多个data frame，然后对每data frame 分别计算，类似于ddply。summarise 可以使用的函数有：min(x), median(x), max(x), quantile(x, p)，计算个数n(), 计算 x 中唯一值的个数n_distinct(), sum(x), mean(x),sum(x > 10), mean(x > 10),sd(x), var(x), iqr(x), mad(x)
group <- group_by(iris_df, Species) # 分组依据
summarise(group, Speciessum = sum(Sepal.Length), Speciesmean=mean(Petal.Length, na.rm = TRUE)) #分组求和## Speciessum Speciesmean
## 1 876 3.8管道操作
%>%与pipeR和magrittr包中%>%操作符一样，用来将上一步产生的对象管道输出为下一步调用的函数的第一个参数。
iris_df %>% group_by(Species) %>%
summarise(total = sum(Sepal.Length)) %>%
arrange(desc(total)) %>%head(5)## total
## 1 876变量查重
通常用select指定需要查重的变量，distinct返回没有重复的数据。
distinct(select(iris_df, Sepal.Length,Species)) #Sepal.Length,Species这两列中没有重复的数据## Source: local data frame [57 x 2]
##
## Sepal.Length Species
## (dbl) (fctr)
## 1 5.1 setosa
## 2 4.9 setosa
## 3 4.7 setosa
## 4 4.6 setosa
## 5 5.0 setosa
## 6 5.4 setosa
## 7 4.4 setosa
## 8 4.8 setosa
## 9 4.3 setosa
## 10 5.8 setosa
## .. ... ...随机抽样
使用sample_n和sample_frac从数据框中随机的返回一些行，sample_n按指定的行数返回，sample_frac按指定的比例返回。
sample_n(iris, 10) #返回10行## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 32 5.4 3.4 1.5 0.4 setosa
## 36 5.0 3.2 1.2 0.2 setosa
## 39 4.4 3.0 1.3 0.2 setosa
## 109 6.7 2.5 5.8 1.8 virginica
## 113 6.8 3.0 5.5 2.1 virginica
## 60 5.2 2.7 3.9 1.4 versicolor
## 92 6.1 3.0 4.6 1.4 versicolor
## 93 5.8 2.6 4.0 1.2 versicolor
## 31 4.8 3.1 1.6 0.2 setosa
## 120 6.0 2.2 5.0 1.5 virginicasample_frac(iris, 0.01) #返回总行数的0.01倍## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 127 6.2 2.8 4.8 1.8 virginica
## 1 5.1 3.5 1.4 0.2 setosa长宽格式数据转换
在wide format中，每一个样本点(subject)自成一行，这一行内记录了这个样本点的所有信息。典型的宽格式数据如下：
data_wide <- read.table(header=TRUE, text='
subject sex control cond1 cond2
1 M 7.9 12.3 10.7
2 F 6.3 10.6 11.1
3 F 9.5 13.1 13.8
4 M 11.5 13.4 12.9
')
data_wide## subject sex control cond1 cond2
## 1 1 M 7.9 12 11
## 2 2 F 6.3 11 11
## 3 3 F 9.5 13 14
## 4 4 M 11.5 13 13long format把wide format中的某几个numerical variables变成了一个factor variable之下的levels，而这几个numerical variables的取值都被集中在了一个变量之下。reshape2包中melt函数把wide format变成long format。
data_long <- melt(data_wide,
id.vars = c('subject', 'sex'),
#ID variables 是指将被保存在long format中的变量,它起到指示样本点的作用
variable.name = 'condition',
value.name = 'measurement')
data_long## subject sex condition measurement
## 1 1 M control 7.9
## 2 2 F control 6.3
## 3 3 F control 9.5
## 4 4 M control 11.5
## 5 1 M cond1 12.3
## 6 2 F cond1 10.6
## 7 3 F cond1 13.1
## 8 4 M cond1 13.4
## 9 1 M cond2 10.7
## 10 2 F cond2 11.1
## 11 3 F cond2 13.8
## 12 4 M cond2 12.9在long format中，每个样本点被拆成了三个行，两个新的变量出现。第一个新变量是一个factor variable，fator levels是wide format中的三个变量。第二个新变量是一个numeric variable，记录的数值对应于wide format中该样本点在control, cond1, cond2三列的取值。除了id.variable之外，其他变量都被变成了long format的形式,数据的长短是相对的，如果把没有转换的sex变量转换掉，数据将变得更长。
melt(data_wide, id.vars = 'subject')## Warning: attributes are not identical across measure variables; they will
## be dropped## subject variable value
## 1 1 sex M
## 2 2 sex F
## 3 3 sex F
## 4 4 sex M
## 5 1 control 7.9
## 6 2 control 6.3
## 7 3 control 9.5
## 8 4 control 11.5
## 9 1 cond1 12.3
## 10 2 cond1 10.6
## 11 3 cond1 13.1
## 12 4 cond1 13.4
## 13 1 cond2 10.7
## 14 2 cond2 11.1
## 15 3 cond2 13.8
## 16 4 cond2 12.9Wide format转换为long format时，最极端的情况是所有变量都转换掉
melt(data_wide, id.vars = NULL)## Warning: attributes are not identical across measure variables; they will
## be dropped## variable value
## 1 subject 1
## 2 subject 2
## 3 subject 3
## 4 subject 4
## 5 sex M
## 6 sex F
## 7 sex F
## 8 sex M
## 9 control 7.9
## 10 control 6.3
## 11 control 9.5
## 12 control 11.5
## 13 cond1 12.3
## 14 cond1 10.6
## 15 cond1 13.1
## 16 cond1 13.4
## 17 cond2 10.7
## 18 cond2 11.1
## 19 cond2 13.8
## 20 cond2 12.9reshape2包中Cast函数把long format变成wide format的函数，dcast针对data.frame，acast针对的是array或matrices。dcast中需要一个formular来说明转换的形式，formular的左边是id variables，右边是一个factor variables，它的factor levels将会在wide format中成为新的variables，这些variables的取值用value.var来指定，最后得到的wide format如下：
data.wide <- dcast(data_long, subject + sex ~ condition, value.var = "measurement")
data.wide## subject sex control cond1 cond2
## 1 1 M 7.9 12 11
## 2 2 F 6.3 11 11
## 3 3 F 9.5 13 14
## 4 4 M 11.5 13 13分类汇总
在比较多组个体或观测时，关注的焦点经常是各组的描述性统计信息，而不是整体的描述性统计信息时，可以使用aggregate()分组获取描述性统计量。 例 epicalc中HW93数据集是1993年泰国南部钩虫感染的调查资料，其中intense变量表示感染的严重程度为有序多分类变量，egp为感染的数量，shoes表示是否穿鞋，agegr是年龄分组,需要计算每个年龄的构虫平均感染钩虫的数量。 使用aggregate()分组获取描述性统计量
data(HW93,package = "epicalc")
aggregate(HW93$epg,by=list(epg=HW93$agegr),mean)## epg mean.HW93$epg
## 1 <15 yrs 1087
## 2 15-59 yrs 908
## 3 60+ yrs 3094注意list(epg=HW93\(age)的使用。如果使用的是list(HW93\)age)，则age列将被标注为Group.1而不是age。如果有多个分组变量，可以使用by=list(name1=groupvar1, name2=groupvar2, … , groupvarN)这样的语句。aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。 doBy包和psych包也提供了分组计算描述性统计量的函数，doBy包中summaryBy()函数的使用格式为：summaryBy(formula,data=dataframe,FUN=function) 其中的formula接受以下的格式：var1+var2+…+varNgrounpvar1+goupvar2+…+groupvarN,在左侧的变量是需要分析的数值型变量，而右侧的变量是类别型的分组变量。function可为任何内建或用户自编的R函数。psych包中的describe.by()函数可计算和describe相同的描述性统计量，只是按照一个或多个分组变量分层，使用psych包中的describe.by()和使用doBy包中的summaryBy()分组计算概述统计量如下，describe.by()函数不允许指定任意函数，所以它的使用范围较窄。若存在一个以上的分组变量，你可以使用list(groupvar1, groupvar2, … , groupvarN)来表示它们。但这仅在分组变量交叉后不出现空白单元时有效。
summaryBy(epg~agegr,data=HW93,FUN=max)## agegr epg.max
## 1 <15 yrs 39123
## 2 15-59 yrs 13223
## 3 60+ yrs 24242describe.by(HW93$epg,HW93$agegr)## Warning: describe.by is deprecated. Please use the describeBy function## group: <15 yrs
## vars n mean sd median trimmed mad min max range skew kurtosis se
## 1 1 259 1087 3279 92 374 136 0 39123 39123 7.3 72 204
## --------------------------------------------------------
## group: 15-59 yrs
## vars n mean sd median trimmed mad min max range skew kurtosis se
## 1 1 331 908 1670 247 506 366 0 13223 13223 3.4 15 92
## --------------------------------------------------------
## group: 60+ yrs
## vars n mean sd median trimmed mad min max range skew kurtosis se
## 1 1 47 3094 6055 690 1504 1023 0 24242 24242 2.6 5.4 883需要使用复杂函数则需要plyr包中的*ply族函数。该函数将这类任务以“分割-应用-结合”这种三步方式进行处理：通过一种或多种factor将数据集进行分割，而后应用某项函数，最后将结果整合回数据集当中。Plyr包中囊括了一整套“ply”函数，其第一个字母表示输入的类型，第二个字母表示输出的类型，输入:array,dataframe,list三种格式，输出: array,dataframe,list,discareded四种格式。plyr包中的ddply()可以得到相同结果。summarize不会提供来自原始数据框中其它列中的任何信息，如果需要列出其它column数据，则可以把“summarize”替换为“transform”，且允许一次应用多个函数。
ddply(.data = HW93,.(agegr),summarize,mean=mean(epg),max=max(epg),min=min(epg))
rate <- function(x){
return(sum(x,na.rm = T)/length(x))
}
ddply(.data = HW93,.(agegr),.fun = function(x){rate(x$epg)})频数表和列联表
table(var1, var2, …, varN) 使用 N 个类别型变量（因子）创建一个 N 维列联表。xtabs(formula, data) 根据一个公式和一个矩阵或数据框创建一个 N 维列联表。prop.table(table, margins) 依margins定义的边际列表将表中条目表示为分数形式。margin.table(table, margins) 依margins定义的边际列表计算表中条目的和addmargins(table, margins) 将概述边margins（默认是求和结果）放入表中。ftable(table) 创建一个紧凑的“平铺”式列联表 ###一维列联表
data(Arthritis,package = "vcd")
pander(head(Arthritis))| ID | Treatment | Sex | Age | Improved |
|---|---|---|---|---|
| 57 | Treated | Male | 27 | Some |
| 46 | Treated | Male | 29 | None |
| 77 | Treated | Male | 30 | None |
| 17 | Treated | Male | 32 | Marked |
| 36 | Treated | Male | 46 | Marked |
| 23 | Treated | Male | 58 | Marked |
mytable<-with(Arthritis,table(Improved))
mytable## Improved
## None Some Marked
## 42 14 28可以用prop.table()将这些频数转化为比例值
prop.table(mytable)## Improved
## None Some Marked
## 0.50 0.17 0.33二维列联表
对于二维列联表，table()函数的使用格式为：table(A,B),其中的A是行变量，B是列变量。xtabs()函数还可使用公式风格的输入创建列联表，格式为：xtabs(~A+B,data=mydata)，其中的mydata是一个矩阵或数据框，要进行交叉分类的变量应出现在公式的右侧（即~符号的右方），以+作为分隔符。若某个变量写在公式的左侧，则其为一个频数向量（在数据已经被表格化时很有用）。
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
mytable## Improved
## Treatment None Some Marked
## Placebo 29 7 7
## Treated 13 7 21可以使用margin.table()和prop.table()函数分别生成边际频数(行和)和比例(行比)。
margin.table(mytable,1)## Treatment
## Placebo Treated
## 43 41prop.table(mytable,1)## Improved
## Treatment None Some Marked
## Placebo 0.67 0.16 0.16
## Treated 0.32 0.17 0.51列和与列比例可以这样计算
margin.table(mytable,2)## Improved
## None Some Marked
## 42 14 28prop.table(mytable,2)## Improved
## Treatment None Some Marked
## Placebo 0.69 0.50 0.25
## Treated 0.31 0.50 0.75各单元格所占比例可用如下语句获取
prop.table(mytable)## Improved
## Treatment None Some Marked
## Placebo 0.345 0.083 0.083
## Treated 0.155 0.083 0.250可以使用addmargins()函数为这些表格添加边际和
addmargins(mytable)## Improved
## Treatment None Some Marked Sum
## Placebo 29 7 7 43
## Treated 13 7 21 41
## Sum 42 14 28 84addmargins(prop.table(mytable))## Improved
## Treatment None Some Marked Sum
## Placebo 0.345 0.083 0.083 0.512
## Treated 0.155 0.083 0.250 0.488
## Sum 0.500 0.167 0.333 1.000在使用addmargins()时，默认是表中所有的变量创建边际和
addmargins(prop.table(mytable,1),2)## Improved
## Treatment None Some Marked Sum
## Placebo 0.67 0.16 0.16 1.00
## Treated 0.32 0.17 0.51 1.00注意 table()函数默认忽略缺失值（NA）。要在频数统计中将NA视为一个有效的类别，请设定参数useNA=“ifany”。
table(Arthritis$Treatment,Arthritis$Improved,useNA = "ifany")##
## None Some Marked
## Placebo 29 7 7
## Treated 13 7 21使用gmodels包中的CrossTable()函数生成二维列联表
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独立性检验；计算期望和（皮尔逊、标准化、调整的标准化）残差；将缺失值作为一种有效值；进行行和列标题的标注;
多维列联表
table()和xtabs()都可以基于三个或更多的类别型变量生成多维列联margin.table()、prop.table()和addmargins()函数也可以推广到多维的情况。另外，ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表
mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis)
mytable## , , 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 5ftable(mytable)## Improved None Some Marked
## Treatment Sex
## Placebo Female 19 7 6
## Male 10 0 1
## Treated Female 6 5 16
## Male 7 2 5margin.table(mytable,c(1,3))#治疗情况（Treatment） × 改善情况（Improved）的边际频数## Improved
## Treatment None Some Marked
## Placebo 29 7 7
## Treated 13 7 21try(detach(package:epicalc))
try(detach(package:mosaic))
try(detach(package:showtext))
try(detach(package:pander))
try(detach(package:PerformanceAnalytics))
try(detach(package:Hmisc))
try(detach(package:pastecs))
try(detach(package:psych))
try(detach(package:plyr))
try(detach(package:doBy))
try(detach(package:vcd))
try(detach(package:gmodels))
try(detach(package:readr))
try(detach(package:stringr))
try(detach(package:dplyr))
try(detach(package:reshape2))