R语言小白学习笔记13—基本统计

笔记链接

学习笔记1—R语言基础.
学习笔记2—高级数据结构.
学习笔记3—R语言读取数据.
学习笔记4—统计图.
学习笔记5—编写R语言函数和简单的控制循环语句.
学习笔记6—分组操作.
学习笔记7—高效的分组操作：dplyr.
学习笔记8—数据迭代.
学习笔记9—数据整理.
学习笔记10—数据重构：Tidyverse.
学习笔记11—字符串操作.
学习笔记12—概率分布.

学习笔记13—基本统计

统计学中最常用的工具有均值、方差、相关系数和t检验，对应R语言的函数有：mean、var、cor和t.test函数。

13.1 概括性统计量

平均值（mean函数和weighted.mean函数）

例：

先生成100个1~100范围内的随机数

这里使用sample函数，其可以从x中均匀取出size大小的数据。

replace=TRUE代表数据可以被抽取多次

> x <- sample(x=1:100, size=100, replace=TRUE)
> x
  [1]  24  38   9  71  29   1  13  93  18  99   9  17  34
 [14]  78  22  17  89  58  67  13  59  34  84   6  89  73
 [27]   1  41  33  21  97  78  15   3  27   4  57  11  31
 [40]  15  87  15   7  85  21   5  66  52  71  55  70  70
 [53]  50  86   9   4  45  90  45  30  62  48  93  75   9
 [66]   3  79  22  94  51  44   9  31  68  24  91  61  87
 [79]  85  30  29  15  55   6  25 100  53  95  64  57  41
 [92]  29  52   3  28  59  51  29  94  33

求均值

> mean(x)
[1] 45.25

在实际情况下可能会遇到缺失值，所以我们随机设置20%的缺失值

> y <- x
> y[sample(x=1:100, size=20, replace=FALSE)] <- NA
> y
  [1]  24  38   9  71  29   1  NA  93  18  99   9  NA  NA
 [14]  78  22  17  89  58  NA  13  59  34  84  NA  89  NA
 [27]   1  41  NA  21  97  78  NA   3  27   4  57  11  31
 [40]  15  87  15   7  85  21   5  66  NA  71  55  NA  NA
 [53]  50  86   9   4  45  90  NA  30  62  48  93  75   9
 [66]   3  NA  NA  94  51  44   9  31  68  24  91  61  87
 [79]  85  30  29  NA  55   6  25 100  53  95  NA  57  NA
 [92]  NA  52  NA  28  59  51  29  NA  33

此时用mean求均值会返回NA，设置参数na.rm=TRUE可以去除数据中缺失值

> mean(y)
[1] NA
> mean(y, na.rm = TRUE)
[1] 46.0375

接下来计算加权平均，需要用weighted.mean函数给出向量的数值和权重

> grades <- c(95, 72, 87, 66)
> weights <- c(1/2, 1/4, 1/8, 1/8)
> mean(grades)
[1] 80
> weighted.mean(x=grades, w=weights)
[1] 84.625

方差（var函数）

例：

> var(x)
[1] 914.7955

标准差（sd函数）

例：

> sd(x)
[1] 30.24559

最大值（max函数），最小值（min），中位数（median）

例：

> min(x)
[1] 1
> max(x)
[1] 100
> median(x)
[1] 42.5
> min(y)
[1] NA
> min(y, na.rm=TRUE)
[1] 1

分位数（quantile函数）

> quantile(x, probs = c(.1, .25, .5, .75, .99))
  10%   25%   50%   75%   99% 
 6.90 17.75 42.50 70.25 99.01 

函数summary可以同时计算均值、最小值、最大值、中位数和部分分位数，且无需设置参数na.rm

> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00   17.75   42.50   45.25   70.25  100.00 

13.2 相关系数和协方差

相关系数：cor函数

这里用ggplot2包中的economics数据

例：

> library(ggplot2)
> head(economics)
# A tibble: 6 x 6
  date         pce    pop psavert uempmed unemploy
  <date>     <dbl>  <dbl>   <dbl>   <dbl>    <dbl>
1 1967-07-01  507. 198712    12.6     4.5     2944
2 1967-08-01  510. 198911    12.6     4.7     2945
3 1967-09-01  516. 199113    11.9     4.6     2958
4 1967-10-01  512. 199311    12.9     4.9     3143
5 1967-11-01  517. 199498    12.8     4.7     3066
6 1967-12-01  525. 199657    11.8     4.8     3018

数据集中pce（个人消费支出）和psavert（个人储蓄率）是反向关系，所以计算其相关系数：

> cor(economics$pce, economics$psavert)
[1] -0.7928546

可以对矩阵用函数cor同时比较多个变量

> cor(economics[, c(2, 4:6)])
                pce    psavert    uempmed   unemploy
pce       1.0000000 -0.7928546  0.7269616  0.6145176
psavert  -0.7928546  1.0000000 -0.3251377 -0.3093769
uempmed   0.7269616 -0.3251377  1.0000000  0.8693097
unemploy  0.6145176 -0.3093769  0.8693097  1.0000000

图形可视化：

这里用GGally包中的函数ggpairs绘制（这个包挺方便的，还是要先安装）：

> GGally::ggpairs(economics[, c(2, 4:6)])

相关关系热图：

#加载reshape2包为了用melt函数
> library(reshape2)
#加载scales包为了绘图
> library(scales)
> econCor <- cor(economics[, c(2, 4:6)])
> econMelt <- melt(econCor, varnames = c("x", "y"), value.name = "Correlation")
> econMelt <- econMelt[order(econMelt$Correlation), ]
> econMelt
          x        y Correlation
2   psavert      pce  -0.7928546
5       pce  psavert  -0.7928546
7   uempmed  psavert  -0.3251377
10  psavert  uempmed  -0.3251377
8  unemploy  psavert  -0.3093769
14  psavert unemploy  -0.3093769
4  unemploy      pce   0.6145176
13      pce unemploy   0.6145176
3   uempmed      pce   0.7269616
9       pce  uempmed   0.7269616
12 unemploy  uempmed   0.8693097
15  uempmed unemploy   0.8693097
1       pce      pce   1.0000000
6   psavert  psavert   1.0000000
11  uempmed  uempmed   1.0000000
16 unemploy unemploy   1.0000000
> ggplot(econMelt, aes(x=x, y=y)) + 
+     geom_tile(aes(fill=Correlation)) +
+     scale_fill_gradient2(low = muted("red"), mid="white",
+                          high="steelblue",
+                          guide=guide_colorbar(ticks=FALSE, barheight = 10),
+                          limits=c(-1, 1)) + 
+     theme_minimal() +
+     labs(x=NULL, y=NULL)

接下来考虑缺失数据的问题

因为cor计算时需要考虑多列数据，所以要用一些特殊的方法处理缺失数据。

例：

先生成数据，数据的第一到三行有缺失值

> m <- c(9, 9, NA, 3, NA, 5, 8, 1, 10, 4)
> n <- c(2, NA, 1, 6, 6, 4, 1, 1, 6, 7)
> p <- c(8, 4, 3, 9, 10, NA, 3, NA, 9, 9)
> q <- c(10, 10, 7, 8, 4, 2, 8, 5, 5, 2)
> r <- c(1, 9, 7, 6, 5, 6, 2, 7, 9, 10)
> theMat <- cbind(m, n, p, q, r)
> theMat
       m  n  p  q  r
 [1,]  9  2  8 10  1
 [2,]  9 NA  4 10  9
 [3,] NA  1  3  7  7
 [4,]  3  6  9  8  6
 [5,] NA  6 10  4  5
 [6,]  5  4 NA  2  6
 [7,]  8  1  3  8  2
 [8,]  1  1 NA  5  7
 [9,] 10  6  9  5  9
[10,]  4  7  9  2 10

方法一：everything

所有列的元素必须不含缺失值，否则返回NA

> cor(theMat, use = "everything")
   m  n  p          q          r
m  1 NA NA         NA         NA
n NA  1 NA         NA         NA
p NA NA  1         NA         NA
q NA NA NA  1.0000000 -0.4242958
r NA NA NA -0.4242958  1.0000000

方法二：all.obs

所有列必须不含缺失值，否则报错

> cor(theMat, use = "all.obs")
Error in cor(theMat, use = "all.obs") : cov/cor中有遗漏值

方法三：complete.obs

仅留下不存在缺失值的行，如果每一行都有缺失值，则会报错

> cor(theMat, use = "complete.obs")
           m          n          p          q          r
m  1.0000000 -0.5228840 -0.2893527  0.2974398 -0.3459470
n -0.5228840  1.0000000  0.8090195 -0.7448453  0.9350718
p -0.2893527  0.8090195  1.0000000 -0.3613720  0.6221470
q  0.2974398 -0.7448453 -0.3613720  1.0000000 -0.9059384
r -0.3459470  0.9350718  0.6221470 -0.9059384  1.0000000

方法四：na.or.complete

仅留下不存在缺失值的行，如果每一行都有缺失值，则会返回NA

> cor(theMat, use = "na.or.complete")
           m          n          p          q          r
m  1.0000000 -0.5228840 -0.2893527  0.2974398 -0.3459470
n -0.5228840  1.0000000  0.8090195 -0.7448453  0.9350718
p -0.2893527  0.8090195  1.0000000 -0.3613720  0.6221470
q  0.2974398 -0.7448453 -0.3613720  1.0000000 -0.9059384
r -0.3459470  0.9350718  0.6221470 -0.9059384  1.0000000

方法五：pairwise.complete

用途较广，它依次比较多对变量，并把两个变量互相之间的缺失值剔除，用余下的数据计算两者相关系数

> cor(theMat, use = "pairwise.complete")
            m           n          p          q
m  1.00000000 -0.02511812 -0.3965859  0.4622943
n -0.02511812  1.00000000  0.8717389 -0.5070416
p -0.39658588  0.87173889  1.0000000 -0.5197292
q  0.46229434 -0.50704163 -0.5197292  1.0000000
r -0.20017222  0.53322585  0.1312506 -0.4242958
           r
m -0.2001722
n  0.5332259
p  0.1312506
q -0.4242958
r  1.0000000

协方差：cov函数

> cov(economics[, c(2, 4:6)])
                  pce      psavert      uempmed
pce      12650851.944 -8359.069071 10618.386190
psavert     -8359.069     8.786360    -3.957847
uempmed     10618.386    -3.957847    16.864531
unemploy  5774578.978 -2422.805358  9431.652268
            unemploy
pce      5774578.978
psavert    -2422.805
uempmed     9431.652
unemploy 6979948.309

13.3 t-检验

t检验数比较两个平均数的差异是否显著

这里用reshape2包中的tips数据集

例：

> data(tips, package = "reshape2")
> head(tips)
  total_bill  tip    sex smoker day   time size
1      16.99 1.01 Female     No Sun Dinner    2
2      10.34 1.66   Male     No Sun Dinner    3
3      21.01 3.50   Male     No Sun Dinner    3
4      23.68 3.31   Male     No Sun Dinner    2
5      24.59 3.61 Female     No Sun Dinner    4
6      25.29 4.71   Male     No Sun Dinner    4

13.3.1 单样本t检验

单样本t检验是检验一个样本平均数与一个已知的总体平均数的差异是否显著

例：

判断变量tip的平均值是否为2.50

> t.test(tips$tip, alternative = "two.sided", mu=2.50)

	One Sample t-test

data:  tips$tip
t = 5.6253, df = 243, p-value = 5.08e-08
alternative hypothesis: true mean is not equal to 2.5
95 percent confidence interval:
 2.823799 3.172758
sample estimates:
mean of x 
 2.998279 

可以看到t检验的步骤，显示了自由度和p值，以及95%的置信区间。

t统计量是被估计均值和假设均值之差与所估计均值标准差的比率。

接下来我们画出数据tip的t分布和t统计量的图

rt函数用来生成标准t分布

> randT <- rt(30000, df=NROW(tips)-1)
> tipTTest <- t.test(tips$tip, alternative = "two.sided", mu=2.50)
> ggplot(data.frame(x=randT)) + 
+     geom_density(aes(x=x), fill="grey", color="grey") +
+     geom_vline(xintercept = tipTTest$statistic) +
+     geom_vline(xintercept = mean(randT) + c(-2, 2)*sd(randT), linetype=2)

13.3.2 两样本t检验

这里我们还用tips中数据比较tip在男女之间是否有差异

例：

首先需要判断其是否服从正态分布

> aggregate(tip ~ sex, data=tips, var)
     sex      tip
1 Female 1.344428
2   Male 2.217424
> shapiro.test(tips$tip)

	Shapiro-Wilk normality test

data:  tips$tip
W = 0.89781, p-value = 8.2e-12

> shapiro.test(tips$tip[tips$sex=="Female"])

	Shapiro-Wilk normality test

data:  tips$tip[tips$sex == "Female"]
W = 0.95678, p-value = 0.005448

> shapiro.test(tips$tip[tips$sex=="Male"])

	Shapiro-Wilk normality test

data:  tips$tip[tips$sex == "Male"]
W = 0.87587, p-value = 3.708e-10

画图：

> ggplot(tips, aes(x=tip, fill=sex)) + 
+     geom_histogram(binwidth = .5, alpha=1/2)

接下来检验两样本方差是否相等

由于检验表明数据不服从正态分布（根据p值），所以不能使用F检验和Bartlett检验。

所以使用非参数Ansari-Bradley检验：

> ansari.test(tip ~ sex, tips)

	Ansari-Bradley test

data:  tip by sex
AB = 5582.5, p-value = 0.376
alternative hypothesis: true ratio of scales is not equal to 1

检验结果表明方差相等，可以用标准的两样本t检验

> t.test(tip ~ sex, data=tips, var.equal=TRUE)

	Two Sample t-test

data:  tip by sex
t = -1.3879, df = 242, p-value = 0.1665
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6197558  0.1074167
sample estimates:
mean in group Female   mean in group Male 
            2.833448             3.089618 

根据检验结果，p值不显著，表明男女工人收到的小费大体相等。

（p值代表的是不接受原假设的最小的显著性水平，可以与选定的显著性水平直接比较。例如取5%的显著性水平，如果P值大于5%，就接受原假设，否则不接受原假设。）

也可以用经验法判断：

> library(plyr)
> tipSummary <- ddply(tips, "sex", summarise,
+                     tip.mean=mean(tip), tip.sd=sd(tip),
+                     Lower=tip.mean - 2*tip.sd/sqrt(NROW(tip)),
+                     Upper=tip.mean + 2*tip.sd/sqrt(NROW(tip)))
> tipSummary
     sex tip.mean   tip.sd    Lower    Upper
1 Female 2.833448 1.159495 2.584827 3.082070
2   Male 3.089618 1.489102 2.851931 3.327304

数据可视化：

> ggplot(tipSummary, aes(x=tip.mean, y=sex)) + geom_point()+
+     geom_errorbarh(aes(xmin=Lower, xmax=Upper), height=.2)

该图显示了根据性别分类的两个样本均值和两倍的标准差区间

13.3.3 两配对样本t检验

R语言中进行配对样本的t检验只需要把t.test中参数paired设置为TRUE即可。

例：

以UsingR包中的父亲和儿子身高数据为例

> data(father.son, package = 'UsingR')
> head(father.son)
   fheight  sheight
1 65.04851 59.77827
2 63.25094 63.21404
3 64.95532 63.34242
4 65.75250 62.79238
5 61.13723 64.28113
6 63.02254 64.24221
> t.test(father.son$fheight, father.son$sheight, paried=TRUE)

	Welch Two Sample t-test

data:  father.son$fheight and father.son$sheight
t = -8.3259, df = 2152.6, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.2317973 -0.7621483
sample estimates:
mean of x mean of y 
 67.68710  68.68407 

检验结果拒绝原假设且认为父亲和儿子有不同的身高

绘制核密度曲线：

> heightDiff <- father.son$fheight - father.son$sheight
> ggplot(father.son, aes(x=fheight-sheight)) +
+     geom_density() +
+     geom_vline(xintercept = mean(heightDiff)) +
+     geom_vline(xintercept = mean(heightDiff) + 
+                    2*c(-1, 1)*sd(heightDiff)/sqrt(nrow(father.son)), linetype=2)

可以看到其分布的均值不为0，且均值的置信区间几乎不包括0

13.4 方差分析

接下来自然是比较多组数据

检验公式：


方差分析用aov函数：

> tipAnova <- aov(tip ~ day - 1, tips)
> tipAnova$coefficients
  dayFri   daySat   daySun  dayThur 
2.734737 2.993103 3.255132 2.771452 

方差分析检验是否有一组与其他组存在差异，但它不能明确哪一组是不同的，所以返回了一个单侧的p值

> summary(tipAnova)
           Df Sum Sq Mean Sq F value Pr(>F)    
day         4 2203.0   550.8   290.1 <2e-16 ***
Residuals 240  455.7     1.9                   
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

由于有一个显著的p值，所以想要查看哪一组与其他组不同，最简单方法为绘制出不同组的估计均值的置信区间图，观察是否有重叠部分。

> tipsByDay <- ddply(tips, "day", plyr::summarize,
+                    tip.mean=mean(tip), tip.sd=sd(tip),
+                    Length=NROW(tip),
+                    tfrac=qt(p=.90, df=Length-1),
+                    Lower=tip.mean - tfrac*tip.sd/sqrt(Length),
+                    Upper=tip.mean + tfrac*tip.sd/sqrt(Length))
> ggplot(tipsByDay, aes(x=tip.mean, y=day)) + geom_point() +
+     geom_errorbarh(aes(xmin=Lower, xmax=Upper), height=.3)