r语言states R语言stat_summary

转载

mob6454cc6e8f43 2024-01-12 17:39:24

文章标签 r语言states r语言descstats 二维缺失值最小值 文章分类 R语言后端开发

文 / 李 博
 描述性统计分析
 频数表和列联表
 相关系数和协方差
 t检验
 非参数统计
7.1描述性统计分析
> myvars
> head(mtcars[myvars])#返回mtcars数据集中的三类
mpg hp wt
Mazda RX4 21.0 110 2.620
Mazda RX4 Wag 21.0 110 2.875
Datsun 710 22.8 93 2.320
Hornet 4 Drive 21.4 110 3.215
Hornet Sportabout 18.7 175 3.440
Valiant 18.1 105 3.460
7.1.1
通过 summary() 计算描述性统计量。
> myvars
> 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
> #包括最小值、最大值、四分位数和数值型变量的均值，
> #以及因子向量和逻辑型向量的频数统计。
通过 sapply() 计算描述性统计量。
sapply(x, FUN, options) 其中的 x 是你的数据框(或矩阵)， FUN 为一个任意的函数。
> options(digits = 3)
> mystats
+ if(na.omit)
+ x
+ m
+ n
+ s
+ skew
+ kurt
+ return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
+ }
> myvars
> sapply(mtcars[myvars],mystats)
mpg hp wt
n 32.000 32.000 32.0000
mean 20.091 146.688 3.2172
stdev 6.027 68.563 0.9785
skew 0.611 0.726 0.4231
kurtosis -0.373 -0.136 -0.0227
7.1.2
通过 Hmisc 包中的 describe() 函数计算描述性统计量
> library(Hmisc)
> myvars
> describe(mtcars[myvars])#变量和观测量、缺失值和唯一值得数目
#平均值、分位数和五个最大、五个最小的值
mtcars[myvars]
3 Variables #变量 32 Observations#观测量
----------------------------------------------------
mpg
n missing distinct Info Mean
32 0 25 0.999 20.09
Gmd .05 .10 .25 .50
6.796 12.00 14.34 15.43 19.20
.75 .90 .95
22.80 30.09 31.30
lowest : 10.4 13.3 14.3 14.7 15.0, highest: 26.0 27.3 30.4 32.4 33.9
----------------------------------------------------
hp
n missing distinct Info Mean
32 0 22 0.997 146.7
Gmd .05 .10 .25 .50
77.04 63.65 66.00 96.50 123.00
.75 .90 .95
180.00 243.50 253.55
lowest : 52 62 65 66 91, highest: 215 230 245 264 335
----------------------------------------------------
wt
n missing distinct Info Mean
32 0 29 0.999 3.217
Gmd .05 .10 .25 .50
1.089 1.736 1.956 2.581 3.325
.75 .90 .95
3.610 4.048 5.293
lowest : 1.513 1.615 1.835 1.935 2.140
highest: 3.845 4.070 5.250 5.345 5.424
----------------------------------------------------
通过 pastecs 包中的 stat.desc() 函数计算描述性统计量。
> options(digits = 3)
> library(pastecs)
> myvars
> stat.desc(mtcars[myvars])
mpg hp wt
nbr.val#所有值 32.00 32.000 32.000
nbr.null#空值 0.00 0.000 0.000
nbr.na#缺失值得数量 0.00 0.000 0.000
min#最小值 10.40 52.000 1.513
max#最大值 33.90 335.000 5.424
range#值域 23.50 283.000 3.911
sum#总和 642.90 4694.000 102.952
median#中位数 19.20 123.000 3.325
mean#平均数 20.09 146.688 3.217
SE.mean#平均数的标准误1.07 12.120 0.173
CI.mean.0.95平均数置信度95%的置信区间
2.17 24.720 0.353
var#方差 36.32 4700.867 0.957
std.dev#标准差 6.03 68.563 0.978
coef.var#变异系数 0.30 0.467 0.304
通过 psych 包中的 describe() 计算描述性统计量
> library(psych)
> myvars
> describe(mtcars[myvars])
vars n mean sd median trimmed mad min
mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40
hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00
wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51
max range skew kurtosis se
mpg 33.90 23.50 0.61 -0.37 1.07
hp 335.00 283.00 0.73 -0.14 12.12
wt 5.42 3.91 0.42 -0.02 0.17
统计量结果包括：非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值得标准误。(忒丰富了，有点变态哈)
7.1.3 分组计算描述性统计量
使用 aggregate() 分组获取描述性统计量
> myvars
> aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)#返回平均值
am mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)#返回标准差
am mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
其中list(am=mtcars$am)的作用是为第一列(am)指定一个标签，若是无指定，则是如下这样：
> myvars
> aggregate(mtcars[myvars],by=list(mtcars$am),mean)
Group.1 mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(mtcars$am),sd)
Group.1 mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
使用 by() 分组计算描述性统计量
by(data, INDICES, FUN),其中data是一个数据框或者矩阵，INDICSE是一个因子或因子组成的列表，定义了分组。FUN是任意函数。
> mystats 
+ if (na.omit)
+ x 
+ m 
+ n 
+ s 
+ skew 
+ kurt 
+ return(c(n=n, mean=m, stdev=s, skew=skew, kurtosis=kurt))
+ }
> dstats
> myvars
> by(mtcars[myvars],mtcars$am,dstats)
mtcars$am: 0
mpg hp wt
n 19.000 19.0000 19.000
mean 17.147 160.2632 3.769
stdev 3.834 53.9082 0.777
skew 0.014 -0.0142 0.976
kurtosis -0.803 -1.2097 0.142
---------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.0000 13.000 13.000
mean 24.3923 126.846 2.411
stdev 6.1665 84.062 0.617
skew 0.0526 1.360 0.210
kurtosis -1.4554 0.563 -1.174
7.1.4 分组计算的扩展
使用 doBy 包中的 summaryBy() 分组计算概述统计量
> library(doBy)
> summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean hp.stdev hp.skew
1 0 19 17.1 3.83 0.0140 -0.803 19 160 53.9 -0.0142
2 1 13 24.4 6.17 0.0526 -1.455 13 127 84.1 1.3599
hp.kurtosis wt.n wt.mean wt.stdev wt.skew wt.kurtosis
1 -1.210 19 3.77 0.777 0.976 0.142
2 0.563 13 2.41 0.617 0.210 -1.174
其中左侧的变量是需要分析的数值型变量(mpg+hp+wt)，右侧(am)是类别型的分组变量。
使用 psych 包中的 describeBy() 分组计算概述统计量
> library(psych)
> myvars
> describeBy(mtcars[myvars],list(am=mtcars$am))
Descriptive statistics by group
am: 0
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80 0.88
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21 12.37
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14 0.18
---------------------------------------------------------------------
am: 1
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46 1.71
hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56 23.31
wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17 0.17
其中describeBy() 函数不允许指定任意函数。
7.2 频数表和列联表(主要是类比型变量，而非描述性统计量)
> library(vcd)
载入需要的程辑包：grid
Warning message:
程辑包‘vcd’是用R版本3.3.3 来建造的
> 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
本节中的数据来自 vcd 包中的 Arthritis 数据集，上例为前几个观测。
7.2.1 生成频数表
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)创建一个紧凑的“平铺”式列联表
1. 一维列联表
> #一维列联表
> mytable
> mytable
Improved
None Some Marked
42 14 28
转化为比例值和百分比：
> prop.table(mytable)
Improved
None Some Marked
0.500 0.167 0.333
> prop.table(mytable)*100
Improved
None Some Marked
50.0 16.7 33.3
2.二维列联表
对于二维列联表，table()函数的使用格式为：mytable
使用mytable 
> mytable
> mytable#对于Arthitis数据集，创建一个自己想要内容的表，然后再利用其它函数进行处理
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
生成边际频数和比例。
>margin.table(mytable,1)#下标1指代table()语句中的第一个变量Treatment
Treatment
Placebo Treated
43 41
> prop.table(mytable,1)# 下标 1 指代 table() 语句中的第一个变量Treatment
Improved
Treatment None Some Marked
Placebo 0.674 0.163 0.163
Treated 0.317 0.171 0.512
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,2)#下标2指的是table中的Improved
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.3452 0.0833 0.0833
Treated 0.1548 0.0833 0.2500
使用 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(prop.table(mytable))
Improved
Treatment None Some Marked Sum
Placebo 0.3452 0.0833 0.0833 0.5119
Treated 0.1548 0.0833 0.2500 0.4881
Sum 0.5000 0.1667 0.3333 1.0000
仅添加各行/各列的和：
> addmargins(prop.table(mytable,1),2)
Improved
Treatment None Some Marked Sum
Placebo 0.674 0.163 0.163 1.000
Treated 0.317 0.171 0.512 1.000
> addmargins(prop.table(mytable,2),1)
Improved
Treatment None Some Marked
Placebo 0.69 0.50 0.25
Treated 0.31 0.50 0.75
Sum 1.00 1.00 1.00
使用 CrossTable 生成二维列联表
> 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 | |
--------------------|-----------|-----------|-----------|-----------|
3.多维列联表
table()和xtabs()都可以基于三个或更多的类别型变量生成多维列联表。
margin.table()、prop.table()和addmargins()函数可以自然地推广到高于二维的情况。
> mytable
> 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 5
边际频数。
> margin.table(mytable,1)#对使用公式 ~Treatement+Sex+Improve 创建的表，
#Treatment需要下标1来引用
Treatment
Placebo Treated
43 41
> margin.table(mytable,2)#Sex需要下标2来引用
Sex
Female Male
59 25
> margin.table(mytable,3)#Improved需要下标3来引用
Improved
None Some Marked
42 14 28
治疗情况( Treatment ) × 改善情况( Improved )的边际频数
> margin.table(mytable,c(1,3))
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
治疗情况( Treatment ) × 性别( Sex )的各类改善情况比例
> ftable(prop.table(mytable,c(1,2)))
Improved None Some Marked
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875
Male 0.9091 0.0000 0.0909
Treated Female 0.2222 0.1852 0.5926
Male 0.5000 0.1429 0.3571
> ftable(addmargins(prop.table(mytable,c(1,2)),3))#以较为紧凑的方式输出多维列联表，
#并添加来边际和。
Improved None Some Marked Sum
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875 1.0000
Male 0.9091 0.0000 0.0909 1.0000
Treated Female 0.2222 0.1852 0.5926 1.0000
Male 0.5000 0.1429 0.3571 1.0000
7.2.2 独立性检验(对生成列联表中的变量是否相关或独立进行检验)
** 1. 卡方独立性检验**(卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度，实际观测值与理论推断值之间的偏离程度就决定卡方值的大小，卡方值越大，越不符合；卡方值越小，偏差越小，越趋于符合，若两个值完全相等时，卡方值就为0，表明理论值完全符合。)
> #卡方独立性检验
> library(vcd)
> mytable
> chisq.test(mytable)#对二维列表的行变量和列变量进行检验
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
> #本例中对治疗情况和改善情况进行检验,P<0.01说明存在一定的相关性。
> mytable
> 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近似算法有可能不准(信息量太少，且有一个小于5的值)
由于P=0.089>0.01，故说明治疗结果和性别之间是不独立的。
2. Fisher精确检验(边界固定的列联表中行和列是相互独立的)
> mytable
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided#替代假设
3. Cochran-Mantel-Haenszel检验(两个名义变量在第三个变量的每一层中都是条件独立的)
> #两个名义变量在第三个变量的每一层中都是条件独立的
> mytable
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2,
p-value = 0.0006647
#表明：分性别来看，用药治疗的患者较接受安慰剂的患者有了更多的改善
7.2.3 相关性的度量(数值越大表明相关性越强)
二维列联表的相关性度量(很是抽象)
> #二维列联表的相关性度量
> library(vcd)
> mytable
> 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
7.3 相关
7.3.1 相关的类型
1. Pearson、Spearman和Kendall相关
使用的函数为：cor(x, use= , method= )，以及cov()函数。
协方差和相关系数
> #协方差和相关系数
> options(digits = 3)
> states
#文盲率、预期寿命、
#谋杀率和高中毕业率数据。
> cov(states)#计算样本协方差(所有变量之间两两计算相关)
Population Income Illiteracy Life Exp Murder HS Grad
Population 19931684 571230 292.868 -407.842 5663.52 -3551.51
Income 571230 377573 -163.702 280.663 -521.89 3076.77
Illiteracy 293 -164 0.372 -0.482 1.58 -3.24
Life Exp -408 281 -0.482 1.802 -3.87 6.31
Murder 5664 -522 1.582 -3.869 13.63 -14.55
HS Grad -3552 3077 -3.235 6.313 -14.55 65.24
> cor(states)# 计算了Pearson积差相关系数
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.0000 0.208 0.108 -0.0681 0.344 -0.0985
Income 0.2082 1.000 -0.437 0.3403 -0.230 0.6199
Illiteracy 0.1076 -0.437 1.000 -0.5885 0.703 -0.6572
Life Exp -0.0681 0.340 -0.588 1.0000 -0.781 0.5822
Murder 0.3436 -0.230 0.703 -0.7808 1.000 -0.4880
HS Grad -0.0985 0.620 -0.657 0.5822 -0.488 1.0000
> cor(states,method = "spearman")# 计算了Spearman等级相关系数
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.000 0.125 0.313 -0.104 0.346 -0.383
Income 0.125 1.000 -0.315 0.324 -0.217 0.510
Illiteracy 0.313 -0.315 1.000 -0.555 0.672 -0.655
Life Exp -0.104 0.324 -0.555 1.000 -0.780 0.524
Murder 0.346 -0.217 0.672 -0.780 1.000 -0.437
HS Grad -0.383 0.510 -0.655 0.524 -0.437 1.000
从表中可以看到收入和高中毕业率之间存在很强的正相关。Pearson相关关系为0.6199，Spearman相关关系为0.510，两者均为极强的正相关，但预期寿命和文盲率存在很强的负相关。Pearson相关关系为-0.588，Spearman相关关系为-0.555。
2. 偏相关(偏相关是指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系)
使用pcor(u, S)函数。
> library(ggm)
> colnames(states)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad"
> pcor(c(1,5,2,3,6),cov(states))#前两个1，5表示要计算相关系数
#的变量下标，其余数值为条件变量的下标。
[1] 0.346#人口和谋杀率之间的相关系数0.346
7.3.2 相关性的显著性检验
使用cor.test(x, y, alternative = , method = )函数，其中的 x 和 y 为要检验相关性的变量， alternative 则用来指定进行双侧检验或单侧检验(取值为 "two.side" 、 "less" 或 "greater" )，而 method 用以指定要计算的相关类型( "pearson" 、"kendall" 或 "spearman" )。
检验某种相关系数的显著性
> #检验某种相关系数的显著性
> cor.test(states[,3],states[,5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 7, df = 50, p-value = 1e-08
alternative hypothesis: true correlation is not equal to 0#替代假设，
95 percent confidence interval:#95%的置信区间
0.528 0.821
sample estimates:#样本估计
cor
0.703
通过 corr.test 计算相关矩阵并进行显著性检验
> library(psych)
> corr.test(states,use = "complete")
Call:corr.test(x = states, use = "complete")
Correlation matrix
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.00 0.21 0.11 -0.07 0.34 -0.10
Income 0.21 1.00 -0.44 0.34 -0.23 0.62
Illiteracy 0.11 -0.44 1.00 -0.59 0.70 -0.66
Life Exp -0.07 0.34 -0.59 1.00 -0.78 0.58
Murder 0.34 -0.23 0.70 -0.78 1.00 -0.49
HS Grad -0.10 0.62 -0.66 0.58 -0.49 1.00
Sample Size
[1] 50
Probability values (Entries above the diagonal are adjusted for multiple tests.)
Population Income Illiteracy Life Exp Murder HS Grad
Population 0.00 0.59 1.00 1.0 0.10 1
Income 0.15 0.00 0.01 0.1 0.54 0
Illiteracy 0.46 0.00 0.00 0.0 0.00 0
Life Exp 0.64 0.02 0.00 0.0 0.00 0
Murder 0.01 0.11 0.00 0.0 0.00 0
HS Grad 0.50 0.00 0.00 0.0 0.00 0
To see confidence intervals of the correlations, print with the short=FALSE option

文 / 李 博
 描述性统计分析
 频数表和列联表
 相关系数和协方差
 t检验
 非参数统计
7.1描述性统计分析
> myvars
> head(mtcars[myvars])#返回mtcars数据集中的三类
mpg hp wt
Mazda RX4 21.0 110 2.620
Mazda RX4 Wag 21.0 110 2.875
Datsun 710 22.8 93 2.320
Hornet 4 Drive 21.4 110 3.215
Hornet Sportabout 18.7 175 3.440
Valiant 18.1 105 3.460
7.1.1
通过 summary() 计算描述性统计量。
> myvars
> 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
> #包括最小值、最大值、四分位数和数值型变量的均值，
> #以及因子向量和逻辑型向量的频数统计。
通过 sapply() 计算描述性统计量。
sapply(x, FUN, options) 其中的 x 是你的数据框(或矩阵)， FUN 为一个任意的函数。
> options(digits = 3)
> mystats
+ if(na.omit)
+ x
+ m
+ n
+ s
+ skew
+ kurt
+ return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
+ }
> myvars
> sapply(mtcars[myvars],mystats)
mpg hp wt
n 32.000 32.000 32.0000
mean 20.091 146.688 3.2172
stdev 6.027 68.563 0.9785
skew 0.611 0.726 0.4231
kurtosis -0.373 -0.136 -0.0227
7.1.2
通过 Hmisc 包中的 describe() 函数计算描述性统计量
> library(Hmisc)
> myvars
> describe(mtcars[myvars])#变量和观测量、缺失值和唯一值得数目
#平均值、分位数和五个最大、五个最小的值
mtcars[myvars]
3 Variables #变量 32 Observations#观测量
----------------------------------------------------
mpg
n missing distinct Info Mean
32 0 25 0.999 20.09
Gmd .05 .10 .25 .50
6.796 12.00 14.34 15.43 19.20
.75 .90 .95
22.80 30.09 31.30
lowest : 10.4 13.3 14.3 14.7 15.0, highest: 26.0 27.3 30.4 32.4 33.9
----------------------------------------------------
hp
n missing distinct Info Mean
32 0 22 0.997 146.7
Gmd .05 .10 .25 .50
77.04 63.65 66.00 96.50 123.00
.75 .90 .95
180.00 243.50 253.55
lowest : 52 62 65 66 91, highest: 215 230 245 264 335
----------------------------------------------------
wt
n missing distinct Info Mean
32 0 29 0.999 3.217
Gmd .05 .10 .25 .50
1.089 1.736 1.956 2.581 3.325
.75 .90 .95
3.610 4.048 5.293
lowest : 1.513 1.615 1.835 1.935 2.140
highest: 3.845 4.070 5.250 5.345 5.424
----------------------------------------------------
通过 pastecs 包中的 stat.desc() 函数计算描述性统计量。
> options(digits = 3)
> library(pastecs)
> myvars
> stat.desc(mtcars[myvars])
mpg hp wt
nbr.val#所有值 32.00 32.000 32.000
nbr.null#空值 0.00 0.000 0.000
nbr.na#缺失值得数量 0.00 0.000 0.000
min#最小值 10.40 52.000 1.513
max#最大值 33.90 335.000 5.424
range#值域 23.50 283.000 3.911
sum#总和 642.90 4694.000 102.952
median#中位数 19.20 123.000 3.325
mean#平均数 20.09 146.688 3.217
SE.mean#平均数的标准误1.07 12.120 0.173
CI.mean.0.95平均数置信度95%的置信区间
2.17 24.720 0.353
var#方差 36.32 4700.867 0.957
std.dev#标准差 6.03 68.563 0.978
coef.var#变异系数 0.30 0.467 0.304
通过 psych 包中的 describe() 计算描述性统计量
> library(psych)
> myvars
> describe(mtcars[myvars])
vars n mean sd median trimmed mad min
mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40
hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00
wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51
max range skew kurtosis se
mpg 33.90 23.50 0.61 -0.37 1.07
hp 335.00 283.00 0.73 -0.14 12.12
wt 5.42 3.91 0.42 -0.02 0.17
统计量结果包括：非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值得标准误。(忒丰富了，有点变态哈)
7.1.3 分组计算描述性统计量
使用 aggregate() 分组获取描述性统计量
> myvars
> aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)#返回平均值
am mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)#返回标准差
am mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
其中list(am=mtcars$am)的作用是为第一列(am)指定一个标签，若是无指定，则是如下这样：
> myvars
> aggregate(mtcars[myvars],by=list(mtcars$am),mean)
Group.1 mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(mtcars$am),sd)
Group.1 mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
使用 by() 分组计算描述性统计量
by(data, INDICES, FUN),其中data是一个数据框或者矩阵，INDICSE是一个因子或因子组成的列表，定义了分组。FUN是任意函数。
> mystats 
+ if (na.omit)
+ x 
+ m 
+ n 
+ s 
+ skew 
+ kurt 
+ return(c(n=n, mean=m, stdev=s, skew=skew, kurtosis=kurt))
+ }
> dstats
> myvars
> by(mtcars[myvars],mtcars$am,dstats)
mtcars$am: 0
mpg hp wt
n 19.000 19.0000 19.000
mean 17.147 160.2632 3.769
stdev 3.834 53.9082 0.777
skew 0.014 -0.0142 0.976
kurtosis -0.803 -1.2097 0.142
---------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.0000 13.000 13.000
mean 24.3923 126.846 2.411
stdev 6.1665 84.062 0.617
skew 0.0526 1.360 0.210
kurtosis -1.4554 0.563 -1.174
7.1.4 分组计算的扩展
使用 doBy 包中的 summaryBy() 分组计算概述统计量
> library(doBy)
> summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean hp.stdev hp.skew
1 0 19 17.1 3.83 0.0140 -0.803 19 160 53.9 -0.0142
2 1 13 24.4 6.17 0.0526 -1.455 13 127 84.1 1.3599
hp.kurtosis wt.n wt.mean wt.stdev wt.skew wt.kurtosis
1 -1.210 19 3.77 0.777 0.976 0.142
2 0.563 13 2.41 0.617 0.210 -1.174
其中左侧的变量是需要分析的数值型变量(mpg+hp+wt)，右侧(am)是类别型的分组变量。
使用 psych 包中的 describeBy() 分组计算概述统计量
> library(psych)
> myvars
> describeBy(mtcars[myvars],list(am=mtcars$am))
Descriptive statistics by group
am: 0
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80 0.88
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21 12.37
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14 0.18
---------------------------------------------------------------------
am: 1
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46 1.71
hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56 23.31
wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17 0.17
其中describeBy() 函数不允许指定任意函数。
7.2 频数表和列联表(主要是类比型变量，而非描述性统计量)
> library(vcd)
载入需要的程辑包：grid
Warning message:
程辑包‘vcd’是用R版本3.3.3 来建造的
> 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
本节中的数据来自 vcd 包中的 Arthritis 数据集，上例为前几个观测。
7.2.1 生成频数表
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)创建一个紧凑的“平铺”式列联表
1. 一维列联表
> #一维列联表
> mytable
> mytable
Improved
None Some Marked
42 14 28
转化为比例值和百分比：
> prop.table(mytable)
Improved
None Some Marked
0.500 0.167 0.333
> prop.table(mytable)*100
Improved
None Some Marked
50.0 16.7 33.3
2.二维列联表
对于二维列联表，table()函数的使用格式为：mytable
使用mytable 
> mytable
> mytable#对于Arthitis数据集，创建一个自己想要内容的表，然后再利用其它函数进行处理
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
生成边际频数和比例。
>margin.table(mytable,1)#下标1指代table()语句中的第一个变量Treatment
Treatment
Placebo Treated
43 41
> prop.table(mytable,1)# 下标 1 指代 table() 语句中的第一个变量Treatment
Improved
Treatment None Some Marked
Placebo 0.674 0.163 0.163
Treated 0.317 0.171 0.512
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,2)#下标2指的是table中的Improved
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.3452 0.0833 0.0833
Treated 0.1548 0.0833 0.2500
使用 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(prop.table(mytable))
Improved
Treatment None Some Marked Sum
Placebo 0.3452 0.0833 0.0833 0.5119
Treated 0.1548 0.0833 0.2500 0.4881
Sum 0.5000 0.1667 0.3333 1.0000
仅添加各行/各列的和：
> addmargins(prop.table(mytable,1),2)
Improved
Treatment None Some Marked Sum
Placebo 0.674 0.163 0.163 1.000
Treated 0.317 0.171 0.512 1.000
> addmargins(prop.table(mytable,2),1)
Improved
Treatment None Some Marked
Placebo 0.69 0.50 0.25
Treated 0.31 0.50 0.75
Sum 1.00 1.00 1.00
使用 CrossTable 生成二维列联表
> 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 | |
--------------------|-----------|-----------|-----------|-----------|
3.多维列联表
table()和xtabs()都可以基于三个或更多的类别型变量生成多维列联表。
margin.table()、prop.table()和addmargins()函数可以自然地推广到高于二维的情况。
> mytable
> 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 5
边际频数。
> margin.table(mytable,1)#对使用公式 ~Treatement+Sex+Improve 创建的表，
#Treatment需要下标1来引用
Treatment
Placebo Treated
43 41
> margin.table(mytable,2)#Sex需要下标2来引用
Sex
Female Male
59 25
> margin.table(mytable,3)#Improved需要下标3来引用
Improved
None Some Marked
42 14 28
治疗情况( Treatment ) × 改善情况( Improved )的边际频数
> margin.table(mytable,c(1,3))
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
治疗情况( Treatment ) × 性别( Sex )的各类改善情况比例
> ftable(prop.table(mytable,c(1,2)))
Improved None Some Marked
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875
Male 0.9091 0.0000 0.0909
Treated Female 0.2222 0.1852 0.5926
Male 0.5000 0.1429 0.3571
> ftable(addmargins(prop.table(mytable,c(1,2)),3))#以较为紧凑的方式输出多维列联表，
#并添加来边际和。
Improved None Some Marked Sum
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875 1.0000
Male 0.9091 0.0000 0.0909 1.0000
Treated Female 0.2222 0.1852 0.5926 1.0000
Male 0.5000 0.1429 0.3571 1.0000
7.2.2 独立性检验(对生成列联表中的变量是否相关或独立进行检验)
** 1. 卡方独立性检验**(卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度，实际观测值与理论推断值之间的偏离程度就决定卡方值的大小，卡方值越大，越不符合；卡方值越小，偏差越小，越趋于符合，若两个值完全相等时，卡方值就为0，表明理论值完全符合。)
> #卡方独立性检验
> library(vcd)
> mytable
> chisq.test(mytable)#对二维列表的行变量和列变量进行检验
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
> #本例中对治疗情况和改善情况进行检验,P<0.01说明存在一定的相关性。
> mytable
> 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近似算法有可能不准(信息量太少，且有一个小于5的值)
由于P=0.089>0.01，故说明治疗结果和性别之间是不独立的。
2. Fisher精确检验(边界固定的列联表中行和列是相互独立的)
> mytable
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided#替代假设
3. Cochran-Mantel-Haenszel检验(两个名义变量在第三个变量的每一层中都是条件独立的)
> #两个名义变量在第三个变量的每一层中都是条件独立的
> mytable
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2,
p-value = 0.0006647
#表明：分性别来看，用药治疗的患者较接受安慰剂的患者有了更多的改善
7.2.3 相关性的度量(数值越大表明相关性越强)
二维列联表的相关性度量(很是抽象)
> #二维列联表的相关性度量
> library(vcd)
> mytable
> 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
7.3 相关
7.3.1 相关的类型
1. Pearson、Spearman和Kendall相关
使用的函数为：cor(x, use= , method= )，以及cov()函数。
协方差和相关系数
> #协方差和相关系数
> options(digits = 3)
> states
#文盲率、预期寿命、
#谋杀率和高中毕业率数据。
> cov(states)#计算样本协方差(所有变量之间两两计算相关)
Population Income Illiteracy Life Exp Murder HS Grad
Population 19931684 571230 292.868 -407.842 5663.52 -3551.51
Income 571230 377573 -163.702 280.663 -521.89 3076.77
Illiteracy 293 -164 0.372 -0.482 1.58 -3.24
Life Exp -408 281 -0.482 1.802 -3.87 6.31
Murder 5664 -522 1.582 -3.869 13.63 -14.55
HS Grad -3552 3077 -3.235 6.313 -14.55 65.24
> cor(states)# 计算了Pearson积差相关系数
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.0000 0.208 0.108 -0.0681 0.344 -0.0985
Income 0.2082 1.000 -0.437 0.3403 -0.230 0.6199
Illiteracy 0.1076 -0.437 1.000 -0.5885 0.703 -0.6572
Life Exp -0.0681 0.340 -0.588 1.0000 -0.781 0.5822
Murder 0.3436 -0.230 0.703 -0.7808 1.000 -0.4880
HS Grad -0.0985 0.620 -0.657 0.5822 -0.488 1.0000
> cor(states,method = "spearman")# 计算了Spearman等级相关系数
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.000 0.125 0.313 -0.104 0.346 -0.383
Income 0.125 1.000 -0.315 0.324 -0.217 0.510
Illiteracy 0.313 -0.315 1.000 -0.555 0.672 -0.655
Life Exp -0.104 0.324 -0.555 1.000 -0.780 0.524
Murder 0.346 -0.217 0.672 -0.780 1.000 -0.437
HS Grad -0.383 0.510 -0.655 0.524 -0.437 1.000
从表中可以看到收入和高中毕业率之间存在很强的正相关。Pearson相关关系为0.6199，Spearman相关关系为0.510，两者均为极强的正相关，但预期寿命和文盲率存在很强的负相关。Pearson相关关系为-0.588，Spearman相关关系为-0.555。
2. 偏相关(偏相关是指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系)
使用pcor(u, S)函数。
> library(ggm)
> colnames(states)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad"
> pcor(c(1,5,2,3,6),cov(states))#前两个1，5表示要计算相关系数
#的变量下标，其余数值为条件变量的下标。
[1] 0.346#人口和谋杀率之间的相关系数0.346
7.3.2 相关性的显著性检验
使用cor.test(x, y, alternative = , method = )函数，其中的 x 和 y 为要检验相关性的变量， alternative 则用来指定进行双侧检验或单侧检验(取值为 "two.side" 、 "less" 或 "greater" )，而 method 用以指定要计算的相关类型( "pearson" 、"kendall" 或 "spearman" )。
检验某种相关系数的显著性
> #检验某种相关系数的显著性
> cor.test(states[,3],states[,5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 7, df = 50, p-value = 1e-08
alternative hypothesis: true correlation is not equal to 0#替代假设，
95 percent confidence interval:#95%的置信区间
0.528 0.821
sample estimates:#样本估计
cor
0.703
通过 corr.test 计算相关矩阵并进行显著性检验
> library(psych)
> corr.test(states,use = "complete")
Call:corr.test(x = states, use = "complete")
Correlation matrix
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.00 0.21 0.11 -0.07 0.34 -0.10
Income 0.21 1.00 -0.44 0.34 -0.23 0.62
Illiteracy 0.11 -0.44 1.00 -0.59 0.70 -0.66
Life Exp -0.07 0.34 -0.59 1.00 -0.78 0.58
Murder 0.34 -0.23 0.70 -0.78 1.00 -0.49
HS Grad -0.10 0.62 -0.66 0.58 -0.49 1.00
Sample Size
[1] 50
Probability values (Entries above the diagonal are adjusted for multiple tests.)
Population Income Illiteracy Life Exp Murder HS Grad
Population 0.00 0.59 1.00 1.0 0.10 1
Income 0.15 0.00 0.01 0.1 0.54 0
Illiteracy 0.46 0.00 0.00 0.0 0.00 0
Life Exp 0.64 0.02 0.00 0.0 0.00 0
Murder 0.01 0.11 0.00 0.0 0.00 0
HS Grad 0.50 0.00 0.00 0.0 0.00 0
To see confidence intervals of the correlations, print with the short=FALSE option

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。