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)