R语言实战之描述性统计分析
下面展示一些 描述性统计分析的R代码语言。
vars <- c("mpg","hp","wt")
head(mtcars[vars])
#创造一个统计的函数列表
#通过sapply()计算描述性统计变量（包括偏度和峰度）
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[vars],mystats)
#sapply(mtcars[myvars],mystats,na.omit=TRUE)单纯忽略缺失值
library(Hmisc)
describe(mtcars[vars])
#通过pastecs包中的stat.desc()函数计算描述性统计量
library(pastecs)
stat.desc(mtcars[vars])
library(psych)
describe(mtcars[vars])
#分组计算描述性统计量
aggregate(mtcars[vars], by = list(am = mtcars$am), mean)
#aggregate()允许在每次调用中使用平均数、标准差这样的单返回值函数。
#使用by()分组计算描述性统计量
dstats <- function(x)sapply(x,mystats)
vars <- c("mpg","hp","wt")
by(mtcars[vars],mtcars$am,dstats)
#使用doBy包中的summaryBy()分组计算概述统计量
library(doBy)
summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats)
#使用psych包分组计算概述统计量
library(psych)
describeBy(mtcars[vars],mtcars$am)
#不允许指定任意函数，所以普适性较低
#psych包中的describeBy()函数可计算和describe()相同的描述性统计量，只按照一个或多个分组的变量分层。
#对比分析
library(doBy)
summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
library(psych)
describeBy(mtcars[vars],list(am=mtcars$am))
#使用reshape包分组计算概述统计量
library(reshape)
dstats1 <- function(x)(c(n = length(x), mean = mean(x), sd = sd(x)))
dfm <- melt(mtcars, measure.vars = vars, id.vars = c("am","cyl"))
cast(dfm, am + cyl + variable ~ . , dstats1)
> vars <- c("mpg","hp","wt")
> head(mtcars[vars])
mpg hp wt
Mazda RX4 21.0 110 2.62
Mazda RX4 Wag 21.0 110 2.88
Datsun 710 22.8 93 2.32
Hornet 4 Drive 21.4 110 3.21
Hornet Sportabout 18.7 175 3.44
Valiant 18.1 105 3.46
>
> #创造一个统计的函数列表
> #通过sapply()计算描述性统计变量（包括偏度和峰度）
> 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[vars],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
> #sapply(mtcars[myvars],mystats,na.omit=TRUE)单纯忽略缺失值
>
> library(Hmisc)
> describe(mtcars[vars]
