title: “使用R进行描述性统计分析”
output:
html_document:
###单组计算描述性统计量
在对数据进行统计分析之前,一般我们需要知道当前数据的描述统计情况,如集中趋势、离散趋势、分布形状。如果数据表中有很多变量,在excel中计算这些统计量的话,要利用公式一个一个进行计算,比较麻烦,在R中可以通过一些简单的函数进行计算。
以车辆路试数据集mtcars为例,这里我们只需要三列数据,每加仑行驶英里数mpg马力hp和车重wt。
data(mtcars)
head(mtcars)
myvars <- c("mpg", "hp","wt")
head(mtcars[myvars])
####方法1:**summary()**函数:
summary()
为R自带的函数,可以返回最大值、最小值、四分位数、平均数6个值。
summary(mtcars[myvars])
####方法2:**apply()或者sapply()函数,
可以计算所选择的任意描述统计量,对于sapply()**函数,其格式为sapply(x, FUN, options),其中x为数据框,FUN为选择的任意函数,例如mean(), sd(), var(), min(), max(), median(), length(), range(), quantile(), fivenum()。得出的结果偏度为+0.61(右偏态),峰度为-0.37,说明较正态分布稍平,
mystats <- function(x, na.omit=FALSE){
if(na.omit)
x <- x[!is.na(x)] #如果x是非空值,则赋值给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, kurt=kurt))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
计算数据的峰度和偏度还可以用e1071包中的**skewness()和kurtosis()**函数来计算。
library(e1071)
mystats <- function(x, na.omit=FALSE){
if(na.omit)
x <- x[!is.na(x)]
skew <- skewness(x)
kurt <- kurtosis(x)
return(c(skew=skew, kurt=kurt))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
####方法3:Hmisc包中的describe()
函数
describe()
相比于方法1,可以返回变量和观测的数目,缺失值、唯一值数目,以及5个最大的值和5个最小的值。
library(Hmisc)
describe(mtcars[myvars])
####方法3:pastecs包中的stat.desc()
函数
基本公式stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
。basic=TRUE
、desc=TRUE
都是默认值,norm=TRUE
(FALSE为默认值)返回正态分布统计量,包括偏度、峰度以及它们的限制性水平和Shapiro-Wilk正态检验结果。该函数还计算了平均数95%的置信区间。
library(pastecs)
stat.desc(mtcars[myvars], norm=TRUE)
####方法4:psych包中的describe()
函数
除了常用的统计量外,可以计算非缺失值的数量、截尾平均数、绝对中位数等。如果截尾平均数和平均数相差太大,说明存在极端值,此时截尾平均数可以更好地反映数据的集中趋势。
注意:hmis包和psych包都含有describe()
函数,但最后载入的程序包优先,hmis包describe()
函数会被掩蔽,如果想要运行先前载入包的函数,可以使用代码hmisc::describe()
,告诉计算机执行哪个函数。
library(psych)
describe(mtcars[myvars])
###分组计算描述性统计量
如果一个变量有多个水平,想计算每个水平数据的描述统计情况,或者有多个被试组,计算每组被试的数据情况,这时可以使用分组计算的方法计算各组数据的描述统计量。
####方法1:**aggregate()**函数
list(am=mtcars$am),am指定了组的列标签。**aggregate()**函数一次只能返回一种统计量
dstats <- function(x)sapply(x, mystats)
myvars <- c("mpg", "hp", "wt")
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)
####方法2:**by()**函数
可以指定任意函数返回多个指定的统计量,其基本公式为by(data, INDICES, FUN)。data为一个数据框;INDICES为一个因子或因子组成的列表,定义分组;FUN为任意函数。这里的dstats调用了前面的mystats函数。
myvars <- c("mpg","hp","wt")
mystats <- function(x, na.omit=FALSE){
if(na.omit)
x <- x[!is.na(x)]
m <- mean(x)
sd <- sd(x)
skew <- skewness(x)
kurt <- kurtosis(x)
return(c(mean=m, SD=sd, skew=skew, kurt=kurt))
}
dstats <- function(x)sapply(x,mystats)
by(mtcars[myvars], mtcars$am, dstats)
####方法3:doBy包中的**summaryBy()**函数
可以指定任意函数返回多个统计量。
library(doBy)
summaryBy(mpg + hp +wt~am, data=mtcars, FUN=mystats)
####方法4:**describeBy()**函数
返回多个统计量,但不可以指定返回哪些统计量
myvars <- c("mpg","hp","wt")
describeBy(mtcars[myvars], list(am=mtcars$am))