前面簡要介紹了R語言的基本數(shù)據(jù)結(jié)構(gòu)和基礎(chǔ)圖形，本節(jié)將簡單介紹如何得到數(shù)據(jù)的描述性統(tǒng)計分析，以及進(jìn)一步了解列聯(lián)表（也即分組）。
1、描述性統(tǒng)計分析
本節(jié)將關(guān)注分析連續(xù)型變量的中心趨勢、變化性和分布形狀的方法。以mtcars數(shù)據(jù)中的mpg、hp和wt為例子對常用的統(tǒng)計方法進(jìn)行說明。
對于數(shù)據(jù)的基礎(chǔ)性統(tǒng)計分析，R中有許多方法，從R自帶的summary()函數(shù)，到三方包中的方法，我們將逐步進(jìn)行介紹。
> df <- mtcars[c('mpg','hp','wt')];head(df)
mpg hp wt
Mazda RX4 21.0 110 2.620
Mazda RX4 Wag 21.0 110 2.875
Datsun 710 22.8 93 2.320
Hornet 4 Drive 21.4 110 3.215
Hornet Sportabout 18.7 175 3.440
Valiant 18.1 105 3.460
> summary(df)
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()函數(shù)提供了最小值、最大值、四分位數(shù)和數(shù)值型變量的均值，以及因子向量和邏輯型向量的頻數(shù)統(tǒng)計。你可以使用apply()函數(shù)或sapply()函數(shù)計算所選擇的任意描述性統(tǒng)計量。函數(shù)fivenum()可返回圖基五數(shù)總括（Tukey’s five-number summary，即最小值、下四分位數(shù)、中位數(shù)、上四分位數(shù)和最大值）。
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,stdenv=s,skew=skew,kurtosis=kurt))
}
> apply(df, 2, mystats)
mpg hp wt
n 32.000000 32.0000000 32.00000000
mean 20.090625 146.6875000 3.21725000
stdenv 6.026948 68.5628685 0.97845744
skew 0.610655 0.7260237 0.42314646
kurtosis -0.372766 -0.1355511 -0.02271075
三方包Hmsic、pastecs和psych中都提供了計算描述性統(tǒng)計的函數(shù)。Hmisc包中的describe()函數(shù)可返回變量和觀測的數(shù)量、缺失值和唯一值的數(shù)目、平均值、分位數(shù)，以及五個最大的值和五個最小的值。
> library(Hmisc)
> describe(df)
df
3 Variables 32 Observations
------------------------------------------------------------------------------------------------------------------------------------------
mpg
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
32 0 25 0.999 20.09 6.796 12.00 14.34 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 .25 .50 .75 .90 .95
32 0 22 0.997 146.7 77.04 63.65 66.00 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 .25 .50 .75 .90 .95
32 0 29 0.999 3.217 1.089 1.736 1.956 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包中有一個名為stat.desc()的函數(shù)，可用來計算種類繁多的描述性統(tǒng)計量。使用方法為stat.desc(x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95),其中的x是一個數(shù)據(jù)框或時間序列。若basic=TRUE（默認(rèn)值），則計算其中所有值、空值、缺失值的數(shù)量，以及最小值、最大值、值域，還有總和。若desc=TRUE（同樣也是默認(rèn)值），則計算中位數(shù)、平均數(shù)、平均數(shù)的標(biāo)準(zhǔn)誤、平均數(shù)置信度為95%的置信區(qū)間、方差、標(biāo)準(zhǔn)差以及變異系數(shù)。最后，若norm=TRUE（不是默認(rèn)的），則返回正態(tài)分布統(tǒng)計量，包括偏度和峰度（以及它們的統(tǒng)計顯著程度）和Shapiro–Wilk正態(tài)檢驗(yàn)結(jié)果。這里使用了p值來計算平均數(shù)的置信區(qū)間（默認(rèn)置信度為0.95）。
> options(digits=2)
> stat.desc(df)
mpg hp wt
nbr.val 32.0 32.00 32.00
nbr.null 0.0 0.00 0.00
nbr.na 0.0 0.00 0.00
min 10.4 52.00 1.51
max 33.9 335.00 5.42
range 23.5 283.00 3.91
sum 642.9 4694.00 102.95
median 19.2 123.00 3.33
mean 20.1 146.69 3.22
SE.mean 1.1 12.12 0.17
CI.mean.0.95 2.2 24.72 0.35
var 36.3 4700.87 0.96
std.dev 6.0 68.56 0.98
coef.var 0.3 0.47 0.30
psych包中也有一個describe()函數(shù)，它可以計算非缺失值的數(shù)量、平均數(shù)、標(biāo)準(zhǔn)差、中位數(shù)、截尾均值、絕對中位差、最小值、最大值、值域、偏度、峰度和平均值的標(biāo)準(zhǔn)誤。
> library(psych)
載入程輯包：‘psych’
The following object is masked from ‘package:Hmisc’:
describe
The following objects are masked from ‘package:ggplot2’:
%+%, alpha
### 這里的describe()函數(shù)被重載了，如果想使用Hmisc中的describe()函數(shù)，可使用Hmisc::describe(df)
> describe(df)
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 32 20.1 6.03 19.2 19.7 5.41 10.4 33.9 23.5 0.61 -0.37 1.07
hp 2 32 146.7 68.56 123.0 141.2 77.10 52.0 335.0 283.0 0.73 -0.14 12.12
wt 3 32 3.2 0.98 3.3 3.1 0.77 1.5 5.4 3.9 0.42 -0.02 0.17
2、分組計算描述性統(tǒng)計量
在前面我們已經(jīng)使用aggregate()函數(shù)進(jìn)行簡單的分組計算，aggregate()的缺點(diǎn)是只能使用平均數(shù)、標(biāo)準(zhǔn)差這樣的單返回值函數(shù)。
若想一次性返回多個統(tǒng)計量，可使用by()函數(shù)，格式為：by(data, Indices, FUN)，其中data是一個數(shù)據(jù)框或矩陣，indices是一個因子或因子組成的列表，定義了分組，F(xiàn)UN則是任意函數(shù)。
### 使用aggregate()函數(shù)分組并獲取單個返回值
> aggregate(df,by=list(mtcars$am),mean)
Group.1 mpg hp wt
1 0 17 160 3.8
2 1 24 127 2.4
###使用by()函數(shù)分組并返回多個統(tǒng)計量
> dstats <- function(x)sapply(x, mystats)
> by(df, mtcars$am, dstats)
mtcars$am: 0
mpg hp wt
n 19.000 19.000 19.00
mean 17.147 160.263 3.77
stdenv 3.834 53.908 0.78
skew 0.014 -0.014 0.98
kurtosis -0.803 -1.210 0.14
-------------------------------------------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.000 13.00 13.00
mean 24.392 126.85 2.41
stdenv 6.167 84.06 0.62
skew 0.053 1.36 0.21
kurtosis -1.455 0.56 -1.17
###這里書中說的dstats()調(diào)用了mystats()函數(shù)，將其應(yīng)用于數(shù)據(jù)框中的每一列中，在通過by()函數(shù)可得到am中每一水平的概括統(tǒng)計量。
###但是個人認(rèn)為應(yīng)該是先分類然后再應(yīng)用統(tǒng)計函數(shù)才對。
doBy包和psych包同樣提供了分祖計算統(tǒng)計量的函數(shù)。doBy包中的summaryBy()函數(shù)使用方法為：summaryBy(formula, data=dataframe, FUN=function)，其中formula接受以下格式：var1+var2+var3+...varN ~ groupvar1+groupvar2+...+groupvarN,在~左側(cè)的變量(varN)是需要分析的數(shù)值型變量，而右側(cè)的變量(groupvarN)是類別型的分組變量。 function可為任何內(nèi)建或用戶自編的R函數(shù)。
> library('doBy')
> summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
am mpg.n mpg.mean mpg.stdenv mpg.skew mpg.kurtosis hp.n hp.mean hp.stdenv hp.skew hp.kurtosis wt.n wt.mean wt.stdenv wt.skew
1 0 19 17 3.8 0.014 -0.8 19 160 54 -0.014 -1.21 19 3.8 0.78 0.98
2 1 13 24 6.2 0.053 -1.5 13 127 84 1.360 0.56 13 2.4 0.62 0.21
wt.kurtosis
1 0.14
2 -1.17
###########
> library(psych)
> describeBy(df)
vars n mean sd median trimmed mad min max range skew kurtosis se
mpg 1 32 20.1 6.03 19.2 19.7 5.41 10.4 33.9 23.5 0.61 -0.37 1.07
hp 2 32 146.7 68.56 123.0 141.2 77.10 52.0 335.0 283.0 0.73 -0.14 12.12
wt 3 32 3.2 0.98 3.3 3.1 0.77 1.5 5.4 3.9 0.42 -0.02 0.17
從上面可以看到doBy中的summaryBy()函數(shù)返回結(jié)果的可讀性比R自帶by()函數(shù)差，但其用起來較為簡單。
psych包中的describeBy()函數(shù)可計算和describe相同的描述性統(tǒng)計量，只是按照一個或多個分組變量分層。但describeBy()函數(shù)不允許指定任意函數(shù)，所以它的普適性較低。若存在一個以上的分組變量，可使用list(groupvar1, groupvar2, ... , groupvarN)來表示它們。
最后，在前面介紹分組方法時，使用過reshape2包，借助reshape2包可以靈活的按組導(dǎo)出描述性統(tǒng)計量（如我前面所言，個人感覺這個包有點(diǎn)難以理解）。但是reshape2相比于reshape有一個很大的缺點(diǎn)，那就是他的聚合函數(shù)只能有一個返回值（雖然reshape2比reshape效率高），就好比aggregate與by。他的使用步驟主要如下：
（1）融合數(shù)據(jù)框: dfm <- melt(dataframe, measure.vars=y, id.vars=g)，其中的dataframe包含著數(shù)據(jù)， y是一個向量，指明了要進(jìn)行概述的數(shù)值型變量（默認(rèn)使用所有變量），而g是由一個或多個分組變量組成的向量；
（2）重鑄數(shù)據(jù)：cast(dfm, groupvar1 + groupvar2 + ... + variable ~ . , FUN), 分組變量以+號分隔，這里的variable只取其字面含義（即僅表示重鑄后數(shù)據(jù)框中的變量variable），而FUN是一個任意函數(shù)。reshape2使用的是dcast()函數(shù)，只能有一個返回值，而cast函數(shù)中卻允許多個返回值。
下面使用mtcars數(shù)據(jù)看看reshape分組聚合的強(qiáng)大之處。
> library(reshape)
載入程輯包：‘reshape’
The following objects are masked from ‘package:reshape2’:
colsplit, melt, recast
### rehape2和reshape包的功能和函數(shù)用法基本相同，除了上面所說的重鑄的統(tǒng)計返回值數(shù)量
> dfm<-melt(mtcars, measure.vars=c('mpg','hp','wt'),id.vars=c('am','cyl'))
> head(dfm)
am cyl variable value
1 1 6 mpg 21.0
2 1 6 mpg 21.0
3 1 4 mpg 22.8
4 0 6 mpg 21.4
5 0 8 mpg 18.7
6 0 6 mpg 18.1
> options(digits=2);cdf<-cast(dfm, am+cyl+variable~., mystats);head(cdf)
am cyl variable n mean stdenv skew kurtosis
1 0 4 mpg 3 22.9 1.45 0.069 -2.3
2 0 4 hp 3 84.7 19.66 -0.380 -2.3
3 0 4 wt 3 2.9 0.41 -0.381 -2.3
4 0 6 mpg 4 19.1 1.63 0.482 -1.9
5 0 6 hp 4 115.2 9.18 -0.094 -2.3
6 0 6 wt 4 3.4 0.12 -0.735 -1.7
3、頻數(shù)表和列聯(lián)表
本部分將著眼于類別型變量的頻數(shù)表和列聯(lián)表。范例中的數(shù)據(jù)來源于vcd包中的Arthritis數(shù)據(jù)集，這份數(shù)據(jù)展示了一項(xiàng)關(guān)于風(fēng)濕性關(guān)節(jié)炎新療法的雙盲臨床試驗(yàn)結(jié)果。
R中提供了多種用于創(chuàng)建頻數(shù)表和列連表的方法，主要常用的函數(shù)如下表：
| 函數(shù) | 描述 |
|---|---|
| table(var1, var2, var3,...,varN) | 使用N個類別型變量（因子）創(chuàng)建一個N維列聯(lián)表 |
| xtabs(formula, data) | 根據(jù)一個公式和一個矩陣或數(shù)據(jù)框創(chuàng)建一個N維列聯(lián)表 |
| pro.table(table, margins) | 依m(xù)argins定義的邊際列表將表中條目表示為分?jǐn)?shù)形式 |
| margin.table(table, margins) | 依m(xù)argins定義的邊際列表計算表中條目的和 |
| addmargins(table,margins) | 將概述邊margins（默認(rèn)是求和結(jié)果）放入表中 |
| ftable(table) | 創(chuàng)建一個緊湊的‘平鋪式’列聯(lián)表 |
接下來，我們將逐個使用以上函數(shù)來探索類別型變量。首先觀察簡單的頻率表，然后是二維列聯(lián)表，最后是多維列聯(lián)表。
3.1 一維列聯(lián)表
一維列聯(lián)表可直接使用table()生成。
> 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 <- table(Arthritis$Improved)
> mytable
None Some Marked
42 14 28
### 產(chǎn)看治療效果的比例
> prop.table(mytable)
None Some Marked
0.50 0.17 0.33
> prop.table(mytable)*100
None Some Marked
50 17 33
###求和
> margin.table(mytable)
[1] 84
###同樣求和，這里將結(jié)果加入了表中
> addmargins(mytable)
None Some Marked Sum
42 14 28 84
> class(mytable)
[1] "table
從上面可以看到，50%的患者都得到了改善。
3.2 二維列聯(lián)表
創(chuàng)建二維列聯(lián)表有兩種方式，使用table(A,B)或xtabs()。xtabs()的使用方式較為復(fù)雜，格式為：xtabs(~A+B, data=mydata)，其中的mydata是一個矩陣或數(shù)據(jù)框?？偟膩碚f，要進(jìn)行交叉分類的變量應(yīng)出現(xiàn)在公式的右側(cè)（即符號的右方），以+作為分隔符。若某個變量寫在公式的左側(cè)，則其為一個頻數(shù)向量（在數(shù)據(jù)已經(jīng)被表格化時很有用），對其進(jìn)行求和運(yùn)算。簡言之，formula公式右邊的應(yīng)該都是類型變量，左邊可選擇加入數(shù)值型變量。
另外值得注意的一點(diǎn)是，table()函數(shù)是默認(rèn)忽略NA的，若需要將NA作為一類，使用useNA='ifany'。
### 使用table()創(chuàng)建二維列聯(lián)表
> table(Arthritis$Treatment,Arthritis$Improved)
None Some Marked
Placebo 29 7 7
Treated 13 7 21
### 使用xtabs()創(chuàng)建二維列聯(lián)表
> mytable <- xtabs(~Treatment+Improved, Arthritis); mytable
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
### 生成邊際頻數(shù)，其中后面的1,2指代xtabs()的formula中變量的下標(biāo)，
### 如 ~Treatment+Improved 中1指的是第一個變量Treatment,2指的是第二個變量Improved
> margin.table(mytable, 1); margin.table(mytable,2)
Treatment
Placebo Treated
43 41
Improved
None Some Marked
42 14 28
### 生成邊際比例
> prop.table(mytable,1)*100;prop.table(mytable,2)*100
Improved
Treatment None Some Marked
Placebo 67 16 16
Treated 32 17 51
Improved
Treatment None Some Marked
Placebo 69 50 25
Treated 31 50 75
### 查看每個單元格所占比例
> prop.table(mytable)
Improved
Treatment None Some Marked
Placebo 0.345 0.083 0.083
Treated 0.155 0.083 0.250
### 添加邊際和
> addmargins(mytable)
Improved
Treatment None Some Marked Sum
Placebo 29 7 7 43
Treated 13 7 21 41
Sum 42 14 28 84
######## 接下來我們再看看，在formula左邊加入一個頻數(shù)變量（這里將Age作為一個頻數(shù)變量）的結(jié)果
> xtabs(Age~Treatment+Improved, Arthritis)
Improved
Treatment None Some Marked
Placebo 1439 402 403
Treated 648 397 1193
### 可以看到關(guān)聯(lián)列表的行列變量沒變，但是其中的數(shù)值變了，這些數(shù)值代表什么意思呢？其實(shí)是每一個分類的和
### 使用aggregate函數(shù)可以簡單實(shí)現(xiàn)相同的結(jié)果
> aggregate(Arthritis$Age,by=list(Treatment=Arthritis$Treatment,Improvement=Arthritis$Improved),sum)
Treatment Improvement x
1 Placebo None 1439
2 Treated None 648
3 Placebo Some 402
4 Treated Some 397
5 Placebo Marked 403
6 Treated Marked 1193
除了R中自帶的table()類函數(shù)外，還可以使用三方包gmodels中的CrossTable()函數(shù)創(chuàng)建二維列聯(lián)表。
> 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 | |
--------------------|-----------|-----------|-----------|-----------|
CrossTable()結(jié)果較為直觀，而且已經(jīng)添加了邊際和和邊際頻率。CrossTable()函數(shù)還有很多選項(xiàng)，可以做許多事情：計算（行、列、單元格）的百分比；指定小數(shù)位數(shù)；進(jìn)行卡方、 Fisher和McNemar獨(dú)立性檢驗(yàn)；計算期望和（皮爾遜、標(biāo)準(zhǔn)化、調(diào)整的標(biāo)準(zhǔn)化）殘差；將缺失值作為一種有效值；進(jìn)行行和列標(biāo)題的標(biāo)注；生成SAS或SPSS風(fēng)格的輸出。
3.3 多維列聯(lián)表
table()和xtabs()都可以基于三個或更多的類別型變量生成多維列聯(lián)。margin.table()、prop.table()和addmargins()函數(shù)可以自然地推廣到高于二維的情況。另外， ftable()函數(shù)可以以一種緊湊而吸引人的方式輸出多維列聯(lián)表。
### 使用table()創(chuàng)建多維列聯(lián)表
> with(Arthritis,table(Treatment, Sex, Improved))
, , Improved = None
Sex
Treatment Female Male
Placebo 19 10
Treated 6 7
, , Improved = Some
Sex
Treatment Female Male
Placebo 7 0
Treated 5 2
, , Improved = Marked
Sex
Treatment Female Male
Placebo 6 1
Treated 16 5
### 使用xtabs創(chuàng)建多維列聯(lián)表
> mytable <- xtabs(~Treatment+Sex+Improved, data=Arthritis)
> mytable
, , Improved = None
Sex
Treatment Female Male
Placebo 19 10
Treated 6 7
, , Improved = Some
Sex
Treatment Female Male
Placebo 7 0
Treated 5 2
, , Improved = Marked
Sex
Treatment Female Male
Placebo 6 1
Treated 16 5
### 使用ftable()對三維列表輸出格式進(jìn)行美化
> ftable(mytable)
Improved None Some Marked
Treatment Sex
Placebo Female 19 7 6
Male 10 0 1
Treated Female 6 5 16
Male 7 2 5
> margin.table(mytable,1)
Treatment
Placebo Treated
43 41
4、小結(jié)
本節(jié)我們介紹了連續(xù)型數(shù)值變量和類別型變量的多種基本統(tǒng)計分析方法，在此對其進(jìn)行一個小結(jié)，推薦一下個人認(rèn)為最好用的方法。
（1）對于連續(xù)型數(shù)值變量的基本統(tǒng)計分析，可利用R內(nèi)置的summary()函數(shù)查看數(shù)據(jù)的描述性統(tǒng)計量（如均值，極值，方差）等，若有特定的要求，可使用apply()、sapply()函數(shù)通過自定義函數(shù)返回想要的結(jié)果。此外三方包Hmisc、pastecs、psych中也都提供了不錯的統(tǒng)計函數(shù)。就輸出結(jié)果的直觀性和功能來說，查看基本統(tǒng)計分析推薦的方法為內(nèi)置summary()函數(shù)和pastecs包中的stat.desc()函數(shù)。
（2）對于分組R同樣提供了多種方法，如果只想對分組結(jié)果查看單個統(tǒng)計分析值（即作用函數(shù)只有單個返回值），aggregate()無疑是最佳選擇，簡單快速。若想對分組結(jié)果進(jìn)行處理返回多個統(tǒng)計值，R同樣有多種方法，這里從輸出直觀性以及使用簡單來說，首推的reshape包（reshape2包無法返回多個統(tǒng)計分析值，因此放棄），其次為doBy包中的summaryBy()函數(shù)和內(nèi)置的by()函數(shù)。
（3）對于類別型變量的分析，一般使用table()和xtabs()即可滿足要求，ftable()能夠是多維關(guān)聯(lián)列表的輸出更為美觀。若想進(jìn)行更高級的分析，則可能需要用到三方包gmodels的CrossTable()函數(shù)。
參考：
R語言實(shí)戰(zhàn).
