R语言学习（六）基本统计分析--上
2019-01-25 本文已影响0人
邱俊辉
整体数据计算描述性统计量
summary()函数
可以使用summary()函数来获取描述性统计量
> 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()函数提供了每个变量的最小值，下四分数，中位数，上四分位数，最大数
sapply()函数
可以使用sapply()函数计算所选的任意描述性统计量，调用格式为
sapply(x,FUN,options)
x是选择的数据框或矩阵，FUN是一个任意的函数
#自定义函数
> 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))
+ }
> 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
在上面的代码中，我们先自编了一个计算观测数，均值，标准差，偏度和峰度的函数，再用sapply()函数，将此函数应用的选择的数据框上
describe()函数
Hmisc包中的describe()函数可以返回变量和观测的数量，缺失值和唯一值的数目，平均值，分位数，以及五个最大的值和五个最小的值
> install.packages("Hmisc")
> library(Hmisc)
> describe(mtcars[myvars])
mtcars[myvars]
3 Variables 32 Observations
-------------------------------------------------------------------------------------
mpg
n missing distinct Info Mean Gmd .05 .10 .25
32 0 25 0.999 20.09 6.796 12.00 14.34 15.43
.50 .75 .90 .95
19.20 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 Gmd .05 .10 .25
32 0 22 0.997 146.7 77.04 63.65 66.00 96.50
.50 .75 .90 .95
123.00 180.00 243.50 253.55
lowest : 52 62 65 66 91, highest: 215 230 245 264 335
-------------------------------------------------------------------------------------
wt
n missing distinct Info Mean Gmd .05 .10 .25
32 0 29 0.999 3.217 1.089 1.736 1.956 2.581
.50 .75 .90 .95
3.325 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
-------------------------------------------------------------------------------------
stac.desc()函数
pastecs包中的stat.desc()函数可以计算数据中的所有值、空值、缺失值的数量、以及最小值、最大值、值域还有总和、中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数
> install.packages("pastecs")
> library(pastecs)
> 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
分组计算描述性统计量
aggregate()函数
> 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
但是aggregate()函数每次只能返回一个统计量要返回多个统计量需要用到by()函数
by()函数
by()函数调用格式为
by(data,INDICES,FUN)
其中data是一个数据框或者矩阵，INDICES是一个因子或因子组成的列表，FUN是函数，如果西药返回多个统计量，需要用户自编函数
> 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))
+ }
> 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
当然你也可以使用别的R包中的计算描述性统计量的函数来自编函数
> by(mtcars[myvars],mtcars$am,describe)#describe是Hmisc函数中的函数，上面可见
mtcars$am: 0
data[x, , drop = FALSE]
3 Variables 19 Observations
-------------------------------------------------------------------------------------
mpg
n missing distinct Info Mean Gmd .05 .10 .25
19 0 16 0.997 17.15 4.454 10.40 12.72 14.95
.50 .75 .90 .95
17.30 19.20 21.76 22.96
Value 10.4 13.3 14.3 14.7 15.2 15.5 16.4 17.3 17.8 18.1 18.7 19.2
Frequency 2 1 1 1 2 1 1 1 1 1 1 2
Proportion 0.105 0.053 0.053 0.053 0.105 0.053 0.053 0.053 0.053 0.053 0.053 0.105
Value 21.4 21.5 22.8 24.4
Frequency 1 1 1 1
Proportion 0.053 0.053 0.053 0.053
-------------------------------------------------------------------------------------
hp
n missing distinct Info Mean Gmd .05 .10 .25
19 0 13 0.993 160.3 63.12 91.7 96.6 116.5
.50 .75 .90 .95
175.0 192.5 233.0 245.0
Value 62 95 97 105 110 123 150 175 180 205 215 230
Frequency 1 1 1 1 1 2 2 2 3 1 1 1
Proportion 0.053 0.053 0.053 0.053 0.053 0.105 0.105 0.105 0.158 0.053 0.053 0.053
Value 245
Frequency 2
Proportion 0.105
-------------------------------------------------------------------------------------
wt
n missing distinct Info Mean Gmd .05 .10 .25
19 0 17 0.996 3.769 0.8082 3.081 3.182 3.438
.50 .75 .90 .95
3.520 3.843 5.269 5.353
Value 2.465 3.150 3.190 3.215 3.435 3.440 3.460 3.520 3.570 3.730 3.780 3.840
Frequency 1 1 1 1 1 3 1 1 1 1 1 1
Proportion 0.053 0.053 0.053 0.053 0.053 0.158 0.053 0.053 0.053 0.053 0.053 0.053
Value 3.845 4.070 5.250 5.345 5.424
Frequency 1 1 1 1 1
Proportion 0.053 0.053 0.053 0.053 0.053
-------------------------------------------------------------------------------------
---------------------------------------------------------------
mtcars$am: 1
data[x, , drop = FALSE]
3 Variables 13 Observations
-------------------------------------------------------------------------------------
mpg
n missing distinct Info Mean Gmd .05 .10 .25
13 0 11 0.995 24.39 7.297 15.48 16.58 21.00
.50 .75 .90 .95
22.80 30.40 32.00 33.00
Value 15.0 15.8 19.7 21.0 21.4 22.8 26.0 27.3 30.4 32.4 33.9
Frequency 1 1 1 2 1 1 1 1 2 1 1
Proportion 0.077 0.077 0.077 0.154 0.077 0.077 0.077 0.077 0.154 0.077 0.077
-------------------------------------------------------------------------------------
hp
n missing distinct Info Mean Gmd .05 .10 .25
13 0 11 0.995 126.8 85.26 59.8 65.2 66.0
.50 .75 .90 .95
109.0 113.0 246.2 292.4
Value 52 65 66 91 93 109 110 113 175 264 335
Frequency 1 1 2 1 1 1 2 1 1 1 1
Proportion 0.077 0.077 0.154 0.077 0.077 0.077 0.154 0.077 0.077 0.077 0.077
-------------------------------------------------------------------------------------
wt
n missing distinct Info Mean Gmd .05 .10 .25
13 0 13 1 2.411 0.7306 1.574 1.659 1.935
.50 .75 .90 .95
2.320 2.780 3.111 3.330
Value 1.513 1.615 1.835 1.935 2.140 2.200 2.320 2.620 2.770 2.780 2.875 3.170
Frequency 1 1 1 1 1 1 1 1 1 1 1 1
Proportion 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077
Value 3.570
Frequency 1
Proportion 0.077
describeBy()函数
describeBy()函数是psych包中的函数，可计算和describe()函数相同的统计量只是按照一个或多个分组变量分层，细心的同学应该会发现describeBy()函数和上面代码中的by(mtcars[myvars],mtcars$am,describe)功效其实一样,只不过美化了一下
调用格式为
describeBy(data,list(a))
其中data为数据框或矩阵，a为分组的向量
> 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
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14
se
mpg 0.88
hp 12.37
wt 0.18
---------------------------------------------------------------
am: 1
vars n mean sd median trimmed mad min max range skew kurtosis
mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46
hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56
wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17
se
mpg 1.71
hp 23.31
wt 0.17
R中还有很多的很多计算描述性统计量的函数，大家可以根据自己爱好选择，这里不做过多介绍了
