「R」基本统计分析
这是来自《R语言实战》的笔记。
因为书中列举的方法和知识点比较多，没必要全都掌握，会一种，其他的了解即可。我就简要地整理一下我觉得重要的吧。
描述性统计分析
R基础包自带summary()函数用于获取描述性统计量，我们调用自带的车辆路试数据集mtcars进行下面相应的展示。
> myvars <- c('mpg', 'hp', 'wt')
> 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
summary函数提供了最小值、最大值、四分位数和数值型变量均值以及因子向量和逻辑型向量的频数统计。
书中提及到了一个叫pastecs的包，包括一个名为stat.desc()的函数，它可以计算种类繁多的描述性统计量。
使用格式：
stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
其中x是一个数据框或时间序列。basic=T计算其中所有值、空值、缺失值的数量，以及最大值、最小值、值域还有总和。desc=T计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。最后，若norm=T（不是默认哈）则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro-Wilk正态检验结果。（是使用p值计算的平均数的置信区间）
> 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
> stat.desc(mtcars[myvars], norm=T)
mpg hp wt
nbr.val 32.0000000 32.00000000 32.00000000
nbr.null 0.0000000 0.00000000 0.00000000
nbr.na 0.0000000 0.00000000 0.00000000
min 10.4000000 52.00000000 1.51300000
max 33.9000000 335.00000000 5.42400000
range 23.5000000 283.00000000 3.91100000
sum 642.9000000 4694.00000000 102.95200000
median 19.2000000 123.00000000 3.32500000
mean 20.0906250 146.68750000 3.21725000
SE.mean 1.0654240 12.12031731 0.17296847
CI.mean.0.95 2.1729465 24.71955013 0.35277153
var 36.3241028 4700.86693548 0.95737897
std.dev 6.0269481 68.56286849 0.97845744
coef.var 0.2999881 0.46740771 0.30412851
skewness 0.6106550 0.72602366 0.42314646
skew.2SE 0.7366922 0.87587259 0.51048252
kurtosis -0.3727660 -0.13555112 -0.02271075
kurt.2SE -0.2302812 -0.08373853 -0.01402987
normtest.W 0.9475647 0.93341934 0.94325772
normtest.p 0.1228814 0.04880824 0.09265499
感觉这个函数统计很全面，基本涵盖了一般统计分析会涉及的基础量和显著性检验结果啊。（用summary加stat.desc是不是够了呢？）
分组计算描述性统计量
可以使用aggregate()函数来分组获取描述性统计量。
> myvars <- c('mpg', 'hp', 'wt')
> aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
am mpg hp wt
1 0 17.14737 160.2632 3.768895
2 1 24.39231 126.8462 2.411000
> aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)
am mpg hp wt
1 0 3.833966 53.90820 0.7774001
2 1 6.166504 84.06232 0.6169816
注意如果使用的是list(am = mtcars$am)这个赋值指定了一个更有帮助的列标签。
因为aggregate()仅允许每次调用中使用单返回值函数，而无法一次返回若干个统计量。为了完成这个任务，可以采用by()函数。格式：
by(data, INDICES, FUN)
其中，data为一个数据框或矩阵，INDICES是一个因子或因子组成的列表，定义了分组，FUN是任意函数。
下面提供一个示例：
> mystats <- function(x, na.omit=FALSE){}
> 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')
> dstats <- function(x) sapply(x, mystats)
> by(mtcars[myvars], mtcars$am, dstats)
mtcars$am: 0
mpg hp wt
n 19.00000000 19.00000000 19.0000000
mean 17.14736842 160.26315789 3.7688947
stdev 3.83396639 53.90819573 0.7774001
skew 0.01395038 -0.01422519 0.9759294
kurtosis -0.80317826 -1.20969733 0.1415676
------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.00000000 13.0000000 13.0000000
mean 24.39230769 126.8461538 2.4110000
stdev 6.16650381 84.0623243 0.6169816
skew 0.05256118 1.3598859 0.2103128
kurtosis -1.45535200 0.5634635 -1.1737358
doBy包和psych包也提供了分组计算描述性统计量的函数，有兴趣自己摸索下（非基础包都要先安装哈）。
频数表和列联表
本节着眼于类别型变量的频数表和列联表，以及相应的独立性检验、相关性的度量、图形化展示结果的方法。除了使用基础安装中的函数，还将使用到vcd包和gmodels包中的函数。
最重要的函数如表：
生成频数表
函数 | 描述 |
|---|---|
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) | 创建一个紧凑的“平铺式”列联表 |
一维列联表
使用table()函数生成简单的频数统计表。
> library(grid)
> library(vcd)
> 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
> mytable <- with(Arthritis, table(Improved))
> mytable
Improved
None Some Marked
42 14 28
# 使用prop.table()将频数转为比例值或百分比
> prop.table(mytable)
Improved
None Some Marked
0.5000000 0.1666667 0.3333333
> prop.table(mytable)*100
Improved
None Some Marked
50.00000 16.66667 33.33333
二维列联表
table()函数使用格式为：
table(A, B) # A为行变量，B为列变量
xtabs()函数可以使用公式风格的输入创建列联表：
mytable <- 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 41
> prop.table(mytable,1) # 1指代table()语句中的第一个变量
Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951
# 使用addmargins()函数为这些表格添加边际和
> addmargins(mytable)
Improved
Treatment None Some Marked Sum
Placebo 29 7 7 43
Treated 13 7 21 41
Sum 42 14 28 84
addmargins(mytable,2)添加行和，addmargins(mytable,1)添加列和。
注意：table()函数默认忽略缺失值（NA）。要在频数统计中将NA视为一个有效的类别，设定参数useNA="ifany"。
使用gmodels包中的CrossTable()函数也可以创建二维列联表，它仿照SAS或SPSS的形式。
> library(gmodels)
> 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 | |
--------------------|-----------|-----------|-----------|-----------|
这个函数有很多功能，参阅help了解详情。
独立性检验
书中描述了3种检验：卡方独立性检验、Fisher精确检验和Cochran-Mantel-Haenszel检验。
卡方独立性检验
使用chisq.test()函数对二联表的行变量和列变量进行卡方独立性检验。
> library(grid)
> library(vcd)
> mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
# 治疗情况和改善情况不独立
> mytable <- xtabs(~ Sex+Improved, data=Arthritis)
> 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近似算法有可能不准
# 性别和改善情况独立
可以通过p-value来看两个变量是否独立。
Fisher精确检验
使用fisher.test()函数进行Fisher精确检验，Fisher检验的原假设是：边界固定的列联表中行和列是相互独立的。
> mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided
Cochran-Mantel-Haenszel检验
使用mantelhaen.test()函数，原假设：两个名义变量在第三个变量中的每一层中都是条件独立的。此检验假设不存在三阶交互作用。
下列代码可以检验治疗情况和改善情况在性别的每一水平下是否独立。(结果表明并不独立)
> mytable <- xtabs(~ Treatment+Improved+Sex, data=Arthritis)
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value =
0.0006647
相关性的度量
上一节的显著性检验评估了是否存在充分的证据以拒绝变量间相互独立的原假设。如果可以拒绝原假设，那么你的兴趣就会自然地转向用以衡量相关性强弱的相关性度量。
vcd包中的assocstats()函数可以用来计算二联表的phi系数、列联系数和Cramer's V系数。
> library(vcd)
> mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
> assocstats(mytable)
X^2 df P(> X^2)
Likelihood Ratio 13.530 2 0.0011536
Pearson 13.055 2 0.0014626
Phi-Coefficient : NA
Contingency Coeff.: 0.367
Cramer's V : 0.394
vcd包提供了一个kappa()函数，可以用于计算混淆矩阵Cohen's kappa值以及加权的kappa值。
vcd包拥有优秀的、用于可视化多维数据集中类别型变量关系的函数，可以绘制马萨克图和关联图。ca包也提供了对应分析函数允许使用多种几何表示。
相关
相关系数可以用来描述定量变量之间的关系。
相关的类型
R可以计算多种相关系数，包括Pearson相关系数、Spearman相关系数、Kendall相关系数、偏相关系数、多分格相关系数和多系列相关系数（具体意义自查）。
Pearson相关系数衡量了两个定量变量之间的线性相关程度。Spearman等级相关系数则衡量分级定序变量之间的相关程度。Kendall's Tau相关系数也是一种非参数的等级相关度量。
cor()计算这三种相关系数，而cov()计算协方差。
> states <- state.x77[, 1:6]
> cov(states)
Population Income Illiteracy Life Exp
Population 19931683.7588 571229.7796 292.8679592 -407.8424612
Income 571229.7796 377573.3061 -163.7020408 280.6631837
Illiteracy 292.8680 -163.7020 0.3715306 -0.4815122
Life Exp -407.8425 280.6632 -0.4815122 1.8020204
Murder 5663.5237 -521.8943 1.5817755 -3.8694804
HS Grad -3551.5096 3076.7690 -3.2354694 6.3126849
Murder HS Grad
Population 5663.523714 -3551.509551
Income -521.894286 3076.768980
Illiteracy 1.581776 -3.235469
Life Exp -3.869480 6.312685
Murder 13.627465 -14.549616
HS Grad -14.549616 65.237894
> cor(states)
Population Income Illiteracy Life Exp Murder
Population 1.00000000 0.2082276 0.1076224 -0.06805195 0.3436428
Income 0.20822756 1.0000000 -0.4370752 0.34025534 -0.2300776
Illiteracy 0.10762237 -0.4370752 1.0000000 -0.58847793 0.7029752
Life Exp -0.06805195 0.3402553 -0.5884779 1.00000000 -0.7808458
Murder 0.34364275 -0.2300776 0.7029752 -0.78084575 1.0000000
HS Grad -0.09848975 0.6199323 -0.6571886 0.58221620 -0.4879710
HS Grad
Population -0.09848975
Income 0.61993232
Illiteracy -0.65718861
Life Exp 0.58221620
Murder -0.48797102
HS Grad 1.00000000
> cor(states, method='spearman')
Population Income Illiteracy Life Exp Murder
Population 1.0000000 0.1246098 0.3130496 -0.1040171 0.3457401
Income 0.1246098 1.0000000 -0.3145948 0.3241050 -0.2174623
Illiteracy 0.3130496 -0.3145948 1.0000000 -0.5553735 0.6723592
Life Exp -0.1040171 0.3241050 -0.5553735 1.0000000 -0.7802406
Murder 0.3457401 -0.2174623 0.6723592 -0.7802406 1.0000000
HS Grad -0.3833649 0.5104809 -0.6545396 0.5239410 -0.4367330
HS Grad
Population -0.3833649
Income 0.5104809
Illiteracy -0.6545396
Life Exp 0.5239410
Murder -0.4367330
HS Grad 1.0000000
cor()函数非常的实用。
偏相关是指控制一个或多个定量变量时，另外两个定量变量之间的相互关系。可以使用ggm包中的pcor()函数来计算，使用前需要安装。
polycor包中的hetcor()函数可以计算一种混合的相关矩阵，有兴趣可以看一下。
相关性的显著性检验
常用的原假设为变量间不相关（总体相关系数为0）。可以使用cor.test()对单个的Pearson、Spearman和Kendall相关系数进行检验。
格式为：cor.test(x, y, alternative = , method =)
> cor.test(states[,3], states[, 5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 6.8479, df = 48, p-value = 1.258e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5279280 0.8207295
sample estimates:
cor
0.7029752
这段代码检验了预期寿命和谋杀率的Pearson相关系数为0的原假设。结果是拒绝原假设（p值很小嘛）。
cor.test()每次只能检验一种相关关系，而psych包中提供的corr.test()可以一次性做更多事情。
我们之前关注了偏相关系数。在多元正态性的假设下，psych包中的pcor.test()可以用来检验控制一个或多个额外变量时两个变量之间的条件独立性。
t检验
这个是统计学领域接触最多的概念了，网上一大堆的解释，我也就不细敲书上的概念了。主要看看怎么用吧。
用的数据集是1960年美国47州的刑罚制度对犯罪率的影响的信息。
独立样本t检验
针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设。这里假设两组数据是独立的，并且从正态总体中抽得。
调用格式：
t.test(y ~ x, data) # 这里y,x都是某个data数据集中的变量
t.test(y1, y2) # 这里就是自己输入的了
这里默认假定方差不相等，并使用Welsh的修正自由度（可以修改）。默认双侧检验（计算都是这样的）。所以注意下参数的设定，没有就直接默认，有特殊要求就help函数的具体参数用法。
> library(MASS)
> t.test(Prob ~ So, data = UScrime)
Welch Two Sample t-test
data: Prob by So
t = -3.8954, df = 24.925, p-value = 0.0006506
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.03852569 -0.01187439
sample estimates:
mean in group 0 mean in group 1
0.03851265 0.06371269
可以依据结果拒绝南方各州和非南方各州拥有相同监禁率的假设（看p值）。
可以在执行t检验之前进行合适的正态化变换（常用就是取log,z变换），不会影响结果。
非独立样本的t检验
假定组间差异呈正态分布。格式：
t.test(y1, y2, paired=TRUE)
> library(MASS)
> with(UScrime, t.test(U1,U2, paired=T))
Paired t-test
data: U1 and U2
t = 32.407, df = 46, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
57.67003 65.30870
sample estimates:
mean of the differences
61.48936
多于两组的情况
使用方差分析（ANOVA），内容单列一章。
组间差异的非参数检验
如果数据无法满足t检验或ANOVA的参数假设，可以转向非参数检验。
两组的比较
若两组数据独立，可以使用Wilcoxon秩和检验（也称为Mann-Whitney U检验）来评估观测是否是从相同的概率分布中抽得的。
调用格式：
wilcox.test(y ~ x, data)
wilcox.test(y1, y2)
默认进行双侧检验，可以添加参数exact来进行精确检验，指定alternative='less'或者greater进行有方向的检验。
Wilcoxon符号秩和检验是非独立样本t检验的一种非参数替代方法。它适用于两组成对数据和无法保证正态性假设的情境。调用格式和上述的Mann-Whitney U检验完全相同，不过还（要）可以添加参数paired = TRUE。
多于两组
如果各组独立，使用Kruskal-Wallis检验：
kruskal.test(y ~ A, data)
如果不独立，使用Friedman检验：
friedman.test(y ~ A | B, data)