使用R进行描述性统计分析(连续性变量)
使用R进行描述性统计分析(连续性变量)
对于描述性统计来说,R可以实现的方法有很多,基础自带的有summary()
函数,还有其他packages,如Hmisc包,pastecs包,psych包提供了计算更多内容的函数。
基础函数
在R中,我们经常使用summary()
函数来计算最大值、最小值、四分位数、均值、频数等等。
data(mtcars)
myvars <- c("mpg", "hp", "wt")
summary(mtcars[myvars])
## --- output------
## NOT RUN:
> 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()
函数就可以得到我们想要的描述性统计量了。不过summary()
函数提供的统计量较少,有时候满足不了我们的需求,那么我们可以使用其他包中提供的函数来进行计算。
其他方法
Hmisc包
Hmisc包是一个包含了很多数据分析函数的包,包括样本量大小的计算,图标绘制,字符串操作,输出为LaTeX及HTML格式的文档等等。在这里可以查看更多详细的信息:
Contains many functions useful for data analysis, high-level graphics, utility operations, functions for computing sample size and power, importing and annotating datasets, imputing missing values, advanced table making, variable clustering, character string manipulation, conversion of R objects to LaTeX and html code, and recoding variables.
在Hmisc包中的describe()
函数提供了数量,缺失值,唯一值的数量,平均数,分位数,**基尼平均值(Geni mean difference, Gmd)**以及五个最大值和最小值:
library(Hmisc)
describe(mtcars[myvars])
## --- output---
## NOT RUN
> describe(mtcars[myvars])
mtcars[myvars] 3 Variables 32 Observations
--------------------------------------------------------------------------------
mpg n missing distinct Info Mean Gmd .05 .10 32 0 25 0.999 20.09 6.796 12.00 14.34 .25 .50 .75 .90 .95 15.43 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 32 0 22 0.997 146.7 77.04 63.65 66.00 .25 .50 .75 .90 .95 96.50 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 32 0 29 0.999 3.217 1.089 1.736 1.956 .25 .50 .75 .90 .95 2.581 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
--------------------------------------------------------------------------------
关于基尼平均值是什么,可以看这里的介绍了解更多相关的内容。
pastecs包
有时候我们想要知道标准差,值域,方差,平均数的95%置信区间,是否符合正态等等结果,那么可能上边的提供的方法无法满足我们的需求。那么可以使用pastecs包中的stat.desc()
函数来计算相关统计量。这个函数将返还一个数据框,这种格式的数据或许比起列表更加容易后续的输出和操作。stat.desc()
函数的主要使用形式是:stat.desc(x,basic = TRUE,desc=TRUE,norm=FALSE,p=0.95)
,其中:
- x:一个数据框对象
- basic:默认为TRUE,计算其中所有值、空值、缺失值的数量,以及最小值、最大值、值域,还有总和
- desc:默认为TRUE,计算中位数、平均数、平均数的标准误、平均数95%置信区间、方差、标准差以及变异系数
- norm:默认为FALSE,计算正态分布统计量,包括偏度和峰度以及它们的统计显著差异和Shapiro-Wilk正态检验结果
library(pastecs)
stat.desc(mtcars[myvars], norm = TRUE, p = 0.95)
## ---output---
## NOT RUN
> stat.desc(mtcars[myvars], norm = TRUE, p = 0.95)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
psych包
在psych包中也提供了一个describe()
函数来计算一般统计量,它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值,绝对中位数、最小值、最大值、值域、偏度、峰度和平均值的标准误:
library(psych)
describe(mtcars[myvars])
## ---output ---
## NOT RUN
> describe(mtcars[myvars])vars n mean sd median trimmed mad min max range skew kurtosis
mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40 33.90 23.50 0.61 -0.37
hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00 335.00 283.00 0.73 -0.14
wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51 5.42 3.91 0.42 -0.02se
mpg 1.07
hp 12.12
wt 0.17
利用sapply()
函数计算描述性统计量
如果对于上述方法提供的结果还不是很满意,那么怎么办?在这种情况下,我们可以考虑使用sapply()
函数来实现我们自定义的统计学描述。关于sapply()
函数,大家是否会回想起以前我们使用的那个lapply()
函数呢?通过查询文档(使用?sapply
)我们可以看到以下的一些信息:
‘lapply’ returns a list of the same length as ‘X’, each element of which is the result of applying ‘FUN’ to the corresponding element of ‘X’.
‘sapply’ is a user-friendly version and wrapper of ‘lapply’ by default returning a vector, matrix or, if ‘simplify = “array”’, an array if appropriate, by applying ‘simplify2array()’. ‘sapply(x, f, simplify = FALSE, USE.NAMES = FALSE)’ is the same as ‘lapply(x, f)’.
这些信息有些难以阅读,但是初步看来,这俩个函数差不多,但是sapply()
函数是一个用户友好版本,而且封装了lapply()
函数,使其返还成向量,矩阵或者数组。不用在意那么多细节,我们看看例子也许就能明白了。
首先,我们需要一个自建的函数来满足我们所需要的统计量:
## -------------------------bulid a func------------------------
mystats <- function(x, na.omit = FALSE){if(na.omit)x <- x[!is.na(x)] # exculde the na datam <- mean(x)n <- length(x)s <- sd(x)skew <- sum((x - m) ^ 3 / s ^ 3) / n # Skewness kurt <- sum((x - m) ^ 4 / s ^ 4) / n # Kurtosisreturn(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt))
}
在这里,我们创建了一个计算均数,数量,标准差,偏度(skewness),峰度(kurtosis)的统计量。接下来我们要对每一个变量进行这些统计量的计算:
sapply(mtcars[myvars], mystats)
## ---output
## NOT RUN
> sapply(mtcars[myvars], mystats)mpg hp wt
n 32.000000 32.0000000 32.0000000
mean 20.090625 146.6875000 3.2172500
stdev 6.026948 68.5628685 0.9784574
skew 0.610655 0.7260237 0.4231465
kurtosis 2.627234 2.8644489 2.9772892
接下来,我们想验证下对于sapply()
函数的理解是否正确,于是打算查看下sapply()
到底返还的是什么类型的对象。想要查看他的结构,首先我们需要把对象进行保存,然后运用str()
函数去查看:
test <- sapply(mtcars[myvars], mystats)
str(test)
## ---output---
## NOT RUN
> str(test)num [1:5, 1:3] 32 20.091 6.027 0.611 2.627 ...- attr(*, "dimnames")=List of 2..$ : chr [1:5] "n" "mean" "stdev" "skew" .....$ : chr [1:3] "mpg" "hp" "wt"
这里我们可以看到他的结构是个二维的5*3的表格,其中所有的数据是num。如果熟悉R语言的数据结构,那么二维的,且每一个元素都相同的表格我们把他定义为矩阵(matrix)。当然我们可以使用is.matrix()
进行验证:
is.matrix(test)
## ---output---
## NOT RUN
> is.matrix(test)
[1] TRUE
由此,我们可以确认对于sapply()
函数的理解大致是准确的:
sapply()
函数是一个用户友好版本,而且封装了lapply()
函数,使其返还成向量,矩阵或者数组。
分组计算
有时候,我们需要的不是计算总体的统计量,而是要计算不同组别的统计量,那么上述的一些方法就不太适用了。我们需要一些其他的方法来实现这个需求。
我们可以使用R自带的aggregate()
函数来计算分组的统计量:
aggregate(mtcars[myvars], by = list(am = mtcars$am), mean)
## ---output---
## NOT RUN
> 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
这里我们分别计算了自动挡(am = 1)组和手动挡(am = 0)组的mpg,hp,wt的均数。aggregate()
函数只能一次计算一个统计量,当需要计算多个统计量的时候需要重复使用,比较麻烦。因此我们需要用其他方式来实现一次多个统计量的计算。
dstatas <- function(x) sapply(x, mystats)
by(mtcars[myvars], mtcars$am, dstatas)
## ---output---
## NOT RUN
> by(mtcars[myvars], mtcars$am, dstatas)
mtcars$am: 0mpg 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 2.19682174 1.79030267 3.1415676
------------------------------------------------------------
mtcars$am: 1mpg hp wt
n 13.00000000 13.000000 13.0000000
mean 24.39230769 126.846154 2.4110000
stdev 6.16650381 84.062324 0.6169816
skew 0.05256118 1.359886 0.2103128
kurtosis 1.54464800 3.563463 1.8262642
在这里,dstatas <- function(x) sapply(x, mystats)
使用了简易的函数写法。使用by()
将数据集分为自动挡和手动挡两组,分别使用函数计算出各个统计量。
除了这种我们自建函数使用by()
函数来进行分组的统计量以外,我们可以使用一些包里提供的方法来计算。
doBy包中的summaryBy()
函数提供了分组计算的功能:
library(doBy)
summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats)
## ---output---
## NOT RUN
> 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
1 0 19 17.14737 3.833966 0.01395038 2.196822 19 160.2632 53.90820
2 1 13 24.39231 6.166504 0.05256118 1.544648 13 126.8462 84.06232hp.skew hp.kurtosis wt.n wt.mean wt.stdev wt.skew wt.kurtosis
1 -0.01422519 1.790303 19 3.768895 0.7774001 0.9759294 3.141568
2 1.35988586 3.563463 13 2.411000 0.6169816 0.2103128 1.826264
psych包中的describeBy()
函数可计算和describe()
相同的描述性统计量,按照一个或多个分组变量进行分层:
library(psych)
describeBy(mtcars[myvars], list(am = mtcars$am))
## ---output---
## NOT RUN
> describeBy(mtcars[myvars], list(am = mtcars$am))Descriptive statistics by group
am: 0vars n mean sd median trimmed mad min max range skew
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98kurtosis se
mpg -0.80 0.88
hp -1.21 12.37
wt 0.14 0.18
------------------------------------------------------------
am: 1vars 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.17se
mpg 1.71
hp 23.31
wt 0.17
describeBy()
函数不允许使用任意指定的函数,所以普适性低,但是胜在于不用自己编写函数,直接就能得出一般的描述性统计量。
输出一般统计描述的表格(连续性变量)
上述的这些方法我们很多都用于数据清洗完的一般性探索中。让我们更加清楚的认清数据的结构,分布等等,为后期的统计建模等等做准备。那么,我们常常见到的医学论文中Table 1的一般统计学描述该如何输出呢?这里我们主要想输出的是连续性变量的平均数和正负标准差(如果符合正态)或者是中位数和四分位数(如果不符合正态)。这里,我们使用自建函数来实现这部分功能:
## 该函数用于一般计数资料的统计学描述
## 当资料符合正态时,使用均数和方差
## 当资料不符合正态时,使用中位数和四分位
library(nortest) # 载入进行正态性检验的包
gl.num.anysis <- function(varnames, source) {# 第一部分:进行正态性检验,得到p值value <- as.vector(as.matrix(source[, varnames]))pvalue <- lillie.test(value)$p.value# 判断p值是否大于0.05,如果大于0.05,为符合正态,并且用‘**’表示符合正态if (pvalue > 0.05) {Mean <- round(mean(value), 4)SD <- round(sd(value), 4)Mean_value <- paste0(Mean,paste0('(',paste(Mean+SD,Mean-SD,sep = '-'),')'),'**')table <- data.frame('Characteristics' = varnames,'Value' = Mean_value)return(table)}# 不符合正态就使用中位数和四分位数,并且用‘*’表示不符合正态else {Median <- paste0(round(median(value),4),paste0('(',paste(round(quantile(value,probs=0.25),4),round(quantile(value,probs=0.75),4),sep = "-"),')'),'*')table <- data.frame('Characteristics' = varnames,'Value' = Median)return(table)}
}
这里,我不对这个自建函数作出过多的解释,主要原理就是首先判断是否符合正态,然后进行分别的运算,使用paste0()
函数进行字符串的操作,黏贴。这里给出下我们示例计算的结果:
gl.num.anysis(varnames = myvars, source = mtcars)
## ---output---
## NOT RUN
> gl.num.anysis(varnames = myvars, source = mtcars)Characteristics Value
1 mpg 19.2(3.69-95.5)*
2 hp 19.2(3.69-95.5)*
3 wt 19.2(3.69-95.5)*
注意,这里存在错误。因为我们在前面使用stat.desc()
函数计算的时候,进行了正态性检验,发现:
> stat.desc(mtcars[myvars], norm = TRUE, p = 0.95)mpg hp wt
normtest.p 0.1228814 0.04880824 0.09265499
可以看到,mpg和wt是符合正态的,为什么在我们上边的结果却判断为了不符合正态,而且每一个值都是19.2。这里的主要原因是函数中的第二步:value <- as.vector(as.matrix(source[, varnames]))
。
这一步将提取出一个变量的数据,并且转换为向量,如果我们使用多个变量,并将其转化为向量,我们会得到这样的结果:
as.vector(as.matrix(mtcars[, myvars]))
## ---output---
## NOT RUN
> as.vector(as.matrix(mtcars[, myvars]))[1] 21.000 21.000 22.800 21.400 18.700 18.100 14.300 24.400 22.800
[10] 19.200 17.800 16.400 17.300 15.200 10.400 10.400 14.700 32.400
[19] 30.400 33.900 21.500 15.500 15.200 13.300 19.200 27.300 26.000
[28] 30.400 15.800 19.700 15.000 21.400 110.000 110.000 93.000 110.000
[37] 175.000 105.000 245.000 62.000 95.000 123.000 123.000 180.000 180.000
[46] 180.000 205.000 215.000 230.000 66.000 52.000 65.000 97.000 150.000
[55] 150.000 245.000 175.000 66.000 91.000 113.000 264.000 175.000 335.000
[64] 109.000 2.620 2.875 2.320 3.215 3.440 3.460 3.570 3.190
[73] 3.150 3.440 3.440 4.070 3.730 3.780 5.250 5.424 5.345
[82] 2.200 1.615 1.835 2.465 3.520 3.435 3.840 3.845 1.935
[91] 2.140 1.513 3.170 2.770 3.570 2.780
所有三个变量全部转换为了一个向量!这显然会得到一个错误的结果。
那么我们应该怎么做?
正确的做法应该是让每一个变量运行一遍这个函数,这里我们将使用我们的老朋友lapply()
函数来解决:
library(plyr)
ldply(lapply(myvars, gl.num.anysis, mtcars))
## ---output---
## NOT RUN
> ldply(lapply(myvars, gl.num.anysis, mtcars))Characteristics Value
1 mpg 20.0906(26.1175-14.0637)**
2 hp 123(96.5-180)*
3 wt 3.2172(4.1957-2.2387)**
这样,我们就得到了正确的结果。
至于想要获得不同分组的统计结果,这里暂时还没有完成,需要大家将数据集切分,同时得到三个表,进行合并,输出,并且在WORD里修改。以后有空在优化吧~
参考文献
- R语言实战,作者:卡巴科弗 ,ISBN: 9787115299901
使用R进行描述性统计分析(连续性变量)相关推荐
- R语言描述性统计分析:使用epiDisplay包的summ函数获取dataframe数据中每个变量的常用统计量、对每个变量进行汇总统计
R语言描述性统计分析:使用epiDisplay包的summ函数获取dataframe数据中每个变量的常用统计量.对每个变量进行汇总统计 目录
- R语言描述性统计分析:相关性分析
R语言描述性统计分析:相关性分析 相关性分析:pearson.spearman.kendall 相关性系数的显著性检验: 偏相关性分析: library(ISwR) attach(thuesen) c ...
- R语言描述性统计分析:假设检验
R语言描述性统计分析:假设检验 单样本t检验: 双样本t检验: 方差齐性检验: 配对样本t检验: 单样本Wilcoxon符号秩检验: 两样本Wilcoxon符号秩检验: daily.intake &l ...
- r 函数返回多个值_第四讲 R描述性统计分析
在"R与生物统计专题"中,我们会从介绍R的基本知识展开到生物统计原理及其在R中的实现.以从浅入深,层层递进的形式在投必得医学公众号更新. 在上一讲中,我们介绍了第三讲 R编程基础- ...
- R语言使用epiDisplay包的tableStack函数基于分组变量生成统计分析表(包含描述性统计分析、假设检验、不同数据使用不同的统计量和假设检验方法)、自定义指定Bartlett检验的p值水平
R语言使用epiDisplay包的tableStack函数基于分组变量生成统计分析表(包含描述性统计分析.假设检验.不同数据使用不同的统计量和假设检验方法).自定义设置assumption.p.val ...
- 在R中对李克特量表(likert)数据进行可视化描述性统计分析,热力图、密度图、柱状图
在R中对李克特量表带数据进行可视化描述性统计分析 李克特量表是一种常用的社会调查问卷模式.常规论文中对多级的李克特量表数据大多计算均值来进行描述性统计分析,但均值较难表现样本整体分布状况,R中like ...
- R语言实战之描述性统计分析
R语言实战之描述性统计分析 下面展示一些 描述性统计分析的R代码语言. vars <- c("mpg","hp","wt") head ...
- R的一些统计分析包工具
[另附]:R语言简明笔记系列 1.1. 描述性统计分析 1.1.1. 描述性统计量的计算 1.1.1.1. summary() > vars<-c("mpg",&quo ...
- python数据分析学习——1.描述性统计分析
描述性分析 描述性统计分析是关于数据的描述和汇总.它使用两种主要方法: 定量方法以数值方式描述和汇总数据. 可视化方法通过图表,曲线图,直方图和其他图形来说明数据. 一般在数据分析的过程中,拿到数据不 ...
- 数据分析师一定要掌握的基础——描述性统计分析
申明:文章内容是作者自己的学习笔记,教学来源是开课吧讲师梁勇老师. 以下博客内容讲解了描述性统计分析的所有知识点,以及利用鸢尾花数据集的分析加强对各个统计量的理解. 数理统计基础-描述性统计分析 1. ...
最新文章
- JavScript中的循环
- 【Hibernate步步为营】--关联映射之多对一
- rust(14)-if let,while let
- Android开发实践:常用NDK命令行参数
- .net mvc html使用方法,C# ASP.NET MVC HtmlHelper用法汇总
- 生命游戏c语言代码,c++生命游戏源码
- 腾讯帝国十八年,被它借鉴过的产品都有哪些?
- 了解Maven的<relativePath/>标签
- PD runner下载和使用教程
- 超越Framer的基础知识
- 利用Gensim训练关于英文维基百科的Word2Vec模型(Training Word2Vec Model on English Wikipedia by Gensim)
- 用lxml的xpath演示爬虫提取笑话集网页其中的标题,url,浏览数,日期,笑话内容
- 常用缓存淘汰策略FIFO、LFU、LRU
- php curl exec 返回值,php curl_exec()函数 CURL获取返回值的方法
- safari下载文件后缀多添加了.exe的解决方法
- 乌班图下载(linux系统)
- 2021年11月中国商品出口总额排行榜:中国贸易顺差717.1亿美元,7个出口最终目的国(地)出口额超过百亿美元(附月榜TOP100详单)
- TypeScript组件化实现弹层播放器
- 使用STM32串口模块配合SerialChart实现虚拟示波器功能
- 2022 最新 Elasticsearch 面试题
热门文章
- 关于申请微信公众号(服务号)的材料和流程
- 大陆、香港、澳门、台湾身份证最全正则校验
- python 爬虫 爬取网易严选全网商品价格评论数据
- mount的挂载远程服务器文件夹
- 360怎样修改wifi服务器地址,360路由器怎么重新设置?
- uni-app背景图片 background-image,支持 base64 格式图片、支持网络路径图片、本地路径背景图片
- c8t6高电平电压_什么是高电平和低电平?
- alex机器人 ser_基于Web Service的机器人远程控制系统设计
- 联想微型计算机拆,联想10064一体机拆机,联想一体机硬盘怎么拆
- e4a php上传,POST上传文件(E4A)