使用R进行描述性统计分析

---

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))