第七章 基本统计分析
我感觉这部分是比较有用的 学会利用函数读取数据，观察数据结构，以便于之后的数据分析工作。
7.1描述性统计分析
这里书里用了R语言内置数据集mtcars来作为演示
7.1.1 summary和sapply
代码：
myvars <- c("mpg","hp","wt","am")
head(mtcars[,myvars])
summary(mtcars[myvars])
结果：
解释：首先定义了一个字符型变量，包含mpg、hp、wt、am四个字符串
用[ ]读取mtcars中名字为这四个的四列，也可以这么写：mtcars[,myvars]是一样的。
summary统计了最值、四分位数、中位数以及平均数。每一列都各自算了
但是只有这些值，要其他更多的描述性统计信息，可以用sapply函数。
格式：sapply(x,FUN,options)
x为数据框，FUN为任意函数，若options被定义了，则会被传递给FUN
代码：
x <- mtcars[,myvars]
sapply(x,FUN = sd)
结果：
也可以sapply fivenum返回图基五数总括（Turkey‘s five-number summary）包括min，下四，中位数，上四，max。
定义option
代码：
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(x,mystats)
注意书写格式 这个就是定义了一个函数，然后将其赋给FUN，sapply来调用这个函数来计算
可以看出来sapply和自建函数结合在一起 可以做很多批量处理的事情！！！
7.1.2 other way
describe函数（Hmics包中）
library(Hmisc)
describe(x)
stat.desc函数
library(pastecs)
stat.desc(x)
计算了这些量
> rownames(stat.desc(x))
[1] "nbr.val" "nbr.null" "nbr.na" "min" "max"
[6] "range" "sum" "median" "mean" "SE.mean"
[11] "CI.mean.0.95" "var" "std.dev" "coef.var"
psych包中的describe函数
library(psych)
describe(x)
计算了这些量
> colnames(describe(x))
[1] "vars" "n" "mean" "sd" "median" "trimmed" "mad"
[8] "min" "max" "range" "skew" "kurtosis" "se"
这里出现了同名的函数，R规则是后来居上，但是Hmisc包中的describe还是存在的，只不过你需要这么读取调用
Hmisc::describe(x)
7.1.3 分组计算描述性统计量
aggregate函数
aggregate(x,by=list(am=x$am),mean)
am是变速箱类型
tips：如果想查看内置数据集每一个变量的含义
?name of the data就行 同理函数也是?function name
> aggregate(x,by=list(am=x$am),mean)
am mpg hp wt am
1 0 17.14737 160.2632 3.768895 0
2 1 24.39231 126.8462 2.411000 1
list(am=x$am)的作用
给分组一个名字
> aggregate(x,by=list(x$am),mean)
Group.1 mpg hp wt am
1 0 17.14737 160.2632 3.768895 0
2 1 24.39231 126.8462 2.411000 1
但是aggregate一次只能返回一个值 mean或者sd，无法一次放回多个统计量
by函数的使用
by(data,INDICES,FUN)
> by(x,INDICES=list(am=x$am),dstats)
am: 0
mpg hp wt am
n 19.00000000 19.00000000 19.0000000 19
mean 17.14736842 160.26315789 3.7688947 0
stdev 3.83396639 53.90819573 0.7774001 0
skew 0.01395038 -0.01422519 0.9759294 NaN
kurtosis -0.80317826 -1.20969733 0.1415676 NaN
-----------------------------------------------------------------
am: 1
mpg hp wt am
n 13.00000000 13.0000000 13.0000000 13
mean 24.39230769 126.8461538 2.4110000 1
stdev 6.16650381 84.0623243 0.6169816 0
skew 0.05256118 1.3598859 0.2103128 NaN
kurtosis -1.45535200 0.5634635 -1.1737358 NaN
7.1.4 其他分组计算的方法
doBY包中的summaryBy函数
> library(doBy)
> summaryBy(mpg+wt+hp~am,data=x,FUN=mystats)
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis wt.n wt.mean wt.stdev
1 0 19 17.14737 3.833966 0.01395038 -0.8031783 19 3.768895 0.7774001
2 1 13 24.39231 6.166504 0.05256118 -1.4553520 13 2.411000 0.6169816
wt.skew wt.kurtosis hp.n hp.mean hp.stdev hp.skew hp.kurtosis
1 0.9759294 0.1415676 19 160.2632 53.90820 -0.01422519 -1.2096973
2 0.2103128 -1.1737358 13 126.8462 84.06232 1.35988586 0.5634635
summaryBy(formula,data=dataframe,FUN=function)
其中formula的格式是 分析变量1+分析变量2+…+分析变量n~分组变量1+…+分组变量n
psych包中的describeBy函数
> library(psych)
> describeBy(x,list(am=x$am))
Descriptive statistics by group
am: 0
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80 0.88
hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21 12.37
wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14 0.18
am 4 19 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN 0.00
-----------------------------------------------------------------
am: 1
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46 1.71
hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56 23.31
wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17 0.17
am 4 13 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00 NaN NaN 0.00
但是describeBy函数没法指定函数 普适性较低
综合下来我觉得summaryBy函数比较简洁好用。
7.2 频数表与列联表
这里书上用了vcd包中的Arthritis数据集作为演示
使用table函数生成一维频数表
> library(vcd)
载入需要的程辑包：grid
> head(Arthritis)
ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked
5 36 Treated Male 46 Marked
6 23 Treated Male 58 Marked
> mytable <- with(Arthritis,table(Improved))
> mytable
Improved
None Some Marked
42 14 28
使用prop.table函数可以变为频率
> prop.table(mytable)
Improved
None Some Marked
0.5000000 0.1666667 0.3333333
> prop.table(mytable)*100
Improved
None Some Marked
50.00000 16.66667 33.33333
二维列联表
使用table函数
> mytable <- with(Arthritis,table(Improved,Treatment))
> mytable
Treatment
Improved Placebo Treated
None 29 13
Some 7 7
Marked 7 21
使用xtabs函数
> mytable <- xtabs(~Treatment+Improved,data = Arthritis)
> mytable
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
用margin.table生成边际频数 然后用prop.table函数生成频数表
> margin.table(mytable,1)
Treatment
Placebo Treated
43 41
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,1)
Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951
> prop.table(mytable,2)
Improved
Treatment None Some Marked
Placebo 0.6904762 0.5000000 0.2500000
Treated 0.3095238 0.5000000 0.7500000
也可以直接用addmargins函数为表格添加边际和行和列
> addmargins(mytable)
Improved
Treatment None Some Marked Sum
Placebo 29 7 7 43
Treated 13 7 21 41
Sum 42 14 28 84
仅添加行或者列的边际和
> addmargins(prop.table(mytable,1),2)
Improved
Treatment None Some Marked Sum
Placebo 0.6744186 0.1627907 0.1627907 1.0000000
Treated 0.3170732 0.1707317 0.5121951 1.0000000
> addmargins(prop.table(mytable,2),1)
Improved
Treatment None Some Marked
Placebo 0.6904762 0.5000000 0.2500000
Treated 0.3095238 0.5000000 0.7500000
Sum 1.0000000 1.0000000 1.0000000
使用gmodels包中的CrossTable函数生成二维列联表
> library(gmodels)
> CrossTable(Arthritis$Treatment,Arthritis$Improved)
Cell Contents
|-------------------------|
| N |
| Chi-square contribution |
| N / Row Total |
| N / Col Total |
| N / Table Total |
|-------------------------|
Total Observations in Table: 84
| Arthritis$Improved
Arthritis$Treatment | None | Some | Marked | Row Total |
--------------------|-----------|-----------|-----------|-----------|
Placebo | 29 | 7 | 7 | 43 |
| 2.616 | 0.004 | 3.752 | |
| 0.674 | 0.163 | 0.163 | 0.512 |
| 0.690 | 0.500 | 0.250 | |
| 0.345 | 0.083 | 0.083 | |
--------------------|-----------|-----------|-----------|-----------|
Treated | 13 | 7 | 21 | 41 |
| 2.744 | 0.004 | 3.935 | |
| 0.317 | 0.171 | 0.512 | 0.488 |
| 0.310 | 0.500 | 0.750 | |
| 0.155 | 0.083 | 0.250 | |
--------------------|-----------|-----------|-----------|-----------|
Column Total | 42 | 14 | 28 | 84 |
| 0.500 | 0.167 | 0.333 | |
--------------------|-----------|-----------|-----------|-----------|
多维列联表不想看了
