R语言统计入门-第五章
第五章 单样本与双样本检验
本章的主题是进行连续型数据的比较,可能是比较两组数据,或者是一组数据与一个预设值。
首先介绍两个函数,用来进行 t 检验的t.test和进行Wilcoxon检验的wilcox.test。它们都能够对单样本、双样本与配对样本进行检验。注意的是,双样本Wilcoxon检验在许多教科书中也叫作Mann-Whitney检验。
5.1 单样本 t 检验
t 检验假设数据来自一个正态分布。在单样本的情况下,假设\(x_1,...,x_n\)来自于服从\(N(\mu,\sigma^2)\)的独立随机变量,其中\(N(\mu,\sigma^2)\)表示均值为\(\mu\)、方差为\(\sigma\)的正态分布。希望对假设\(\mu=\mu_0\)进行检验。我们能用样本的均值\(\overline{x}\)与标准差\(s\)来估计真实参数\(\mu\)与\(\sigma\)。
下面讲下均值的标准误,即SEM(Standard Error of the Mean),对\(n\)个均值为\(\mu\)、方差为\(\sigma\)的随机变量求平均值,SEM便用来描述这个平均值的波动性,它的表达式是: \[ SEM=\sigma/\sqrt{n} \] 这个值能告诉我们观察到的均值偏离真实值多远是比较合适的。对于服从正态分布的数据,有一条一般性准则:有95%的数据会落在\(\mu\pm2\sigma\)这个区间里。所以如果\(\mu_0\)是真实的平均值,那么\(\overline{x}\)就应该落在\(\mu_0\pm2SEM\)中。也就是说,可以通过计算: \[ t=\frac{\overline{x}-\mu_0}{SEM} \] 来判断 t 是否落在了一个接受域中。t 应该以一定的概率落在这个接受域之外,这个概率被称为显著性水平。显著性水平一般被设为5%,此时接受域大致是-2到2的区间。
如果 t 落在接受域之外,那么就在预设的显著性水平上拒绝零假设。另一种(等价的)的方法是计算 p 值,它指的是得到一个绝对值上大于或等于当前 t 值的概率,我们能在 p 值小于显著性水平的情况下拒绝零假设。
下面是一个具体例子,反映11个女性的每日摄入能量(千焦)记录。
daily.intake <- c(5260,5470,6180,6390,6515,6805,7515,7515,8230,8770)
#先看下简单的描述性统计
mean(daily.intake)
## [1] 6865
sd(daily.intake)
## [1] 1139.213
quantile(daily.intake)
## 0% 25% 50% 75% 100%
## 5260.0 6232.5 6660.0 7515.0 8770.0
#现在想检验一下这些数据是否与推荐值7725千焦相差甚远,假设这些数据来自于正态分布
#那么目的就是检验这个分布是否满足μ=7725
#下面用t.test进行
t.test(daily.intake,mu=7725)
##
## One Sample t-test
##
## data: daily.intake
## t = -2.3872, df = 9, p-value = 0.04074
## alternative hypothesis: true mean is not equal to 7725
## 95 percent confidence interval:
## 6050.056 7679.944
## sample estimates:
## mean of x
## 6865
#t.test函数还有其他参数。这里用mu来指定零假设的均值(默认mu=0)
#此外,还可以用alternative="greater"来检验均值是否大于μ
#或者是alternative="less"来检验均值是否小于μ
#还可以通过conf.level=0.99来得到一个99%置信区间
下面是对这个结果的解释:
One Sample t-test
这是对所做检验的描述。注意看函数的调用格式,t.test会自动发现用户需要进行一个单样本检验
data:daily.intake
这个是被检验的数据。
t=-2.8208, df=10, p-value=0.01814
这里显示了 t 统计量,相应的自由度df以及准确的 p 值。这里不用拿着一张 t 分布表来查这个统计量落在了那两个分位数之间,就能够立刻看到 p <0.05,于是在自定义的5%显著性水平下,数据显著偏离了原假设中的均值7725千焦。
alternative hypothesis:true mean is not equal to 7725
这里有两个重要的信息:(a)原假设中均值;(b)这是一个双边检验(“not equal to 7725”)。
95 percent confidence interval:
5986.348 7520.925
这是真实均值的95%置信区间:这个区间是一个集合,如果原假设的值来自于这个集合,那么数据便不会显著的偏离。对 t 检验的步骤进行逆向求解,可以得到使得 t 统计量会落在接受域中的一组\(\mu_0\)。 \[ \overline{x}-t_{0.975}(f)\times SEM\lt\mu\lt\overline{x}+t_{0.975}(f)\times SEM \] sample estimates:
mean of x
6753.636
位于最后的这部分是观测值的均值,也就是对真实均数的(点)估计。
5.2 Wilcoxon符号秩检验
t 检验在用于那些不来自于正态分布的数据时比较稳定,尤其是在大样本情况下。不过有的时候会避免做出这种假设(即非正态分布数据),这时最好用不依赖于分布的方法,这些方法通常都把数据替换成了相应的顺序统计量。
对Wilcoxon秩和检验的实际应用基本上与 t 检验相同。
#与t.test一样,wilcox.test有mu和alternative两个参数
#此外,还有correct这个参数,用来指示是否需要进行连续型修正,比如correct=F表示不用
#还有exact,用来指示是否需要进行精确的计算,也用TRUE和FALSE这两个逻辑值来控制
wilcox.test(daily.intake, mu=7725)
## Warning in wilcox.test.default(daily.intake, mu = 7725): cannot compute
## exact p-value with ties
##
## Wilcoxon signed rank test with continuity correction
##
## data: daily.intake
## V = 8, p-value = 0.05263
## alternative hypothesis: true location is not equal to 7725
这里比t.test的输出更短些,因为一个分参数检验不会出现类似于参数估计以及置信边界等概念。
结果里面的V是检验统计量,表示的是正数对应的秩和,在这个例子中,因为存在两个7515,所以p值是通过正态近似的方法来计算的。
5.3 两样本 t 检验
两样本 t 检验主要用来检验两个样本是否来自于均值相等的分布。
两样本检验和单样本检验的理论基础差别不大。现在从两个组别中抽出数据\(x_{11}...,x_{1n_1}\)和\(x_{21},...,x_{2n_2}\),假设它们是从\(N(\mu_1,\sigma_1^2)\)与\(N(\mu_2,\sigma_2^2)\)两个分布中抽取的样本,并希望检验零假设\(\mu_1=\mu_2\)。接着计算: \[ t=\frac{\overline{x}_2-\overline{x}_1}{SEDM} \] SEDM(Standard Error of Difference of Means)是均值差的标准误,被定义为: \[ SEDM=\sqrt{SEM_1^2+SEM^2_2} \] 对于是否假设两组数据的方差相等,SEDM有两种相应计算的方法。“经典”的方法假设方差相等。在这个方法下,先根据两组的标准差计算一个综合性的 s ,然后将它代入SEM。零假设下的 t 值服从自由度为\(n_1+n_2-2\)的 t 分布。
另一种方法由Welch提出,即以两组各自的标准差\(s_1\)和\(s_2\)来计算SEM。这个方法算出来 t 统计量已经不满足 t 分布了,不过可以通过一个 t 分布来近似。这个近似分布的自由度能够通过\(s_1,s_2\)与样本量来得到,一般不是整数。
回头看看每日消耗鞥能量的数据(见1.2.14小节),来比较一下两组女性每日消耗的能量是否有差别。
library(ISwR)
attach(energy)
energy
## expend stature
## 1 9.21 obese
## 2 7.53 lean
## 3 7.48 lean
## 4 8.08 lean
## 5 8.09 lean
## 6 10.15 lean
## 7 8.40 lean
## 8 10.88 lean
## 9 6.13 lean
## 10 7.90 lean
## 11 11.51 obese
## 12 12.79 obese
## 13 7.05 lean
## 14 11.85 obese
## 15 9.97 obese
## 16 7.48 lean
## 17 8.79 obese
## 18 9.69 obese
## 19 9.68 obese
## 20 7.58 lean
## 21 9.19 obese
## 22 8.11 lean
这个数据框的两列包含了所有我们需要的信息。属性变量satrure包含了分组信息,而数值变量expend包含了以兆焦耳为单位的能量消耗。我们只要传递一个模型方程,就能通过R中的t.test和wilcox.test来分析这种格式的数据。
我们的目的是两组的水平是否有差异,所以我们用一个 t 检验。
#~波浪号是致命expend是通过stature来描述的
t.test(expend~stature)
##
## Welch Two Sample t-test
##
## data: expend by stature
## t = -3.8555, df = 15.919, p-value = 0.001411
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.459167 -1.004081
## sample estimates:
## mean in group lean mean in group obese
## 8.066154 10.297778
输出的内容和单样本检验中的基本一致。其中均值之差的置信区间不包含0,与 p 值的结果相统一,意味着在5%的置信水平下差异是显著的。
R中的 t 检验默认采用Welch的变种算法,当不假设两组方差相等时就应该用它,这也会导致非整数的自由度。
为进行平常的 t 检验,应该明确方差相等这个假设,可以通过使参数var.equal=T来达到这一点:
t.test(expend~stature, var.equal=T)
##
## Two Sample t-test
##
## data: expend by stature
## t = -3.9456, df = 20, p-value = 0.000799
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.411451 -1.051796
## sample estimates:
## mean in group lean mean in group obese
## 8.066154 10.297778
5.4 比较方差
虽然在R中不需要假设方差相等也能进行两样本 t 检验,但仍可能对这个假设本身是否正确感兴趣。为此,R提供了var.test函数来做到这一点。这个函数的主要功能是对两样本的方差进行 F 检验。它的使用方法和t.test一样:
var.test(expend~stature)
##
## F test to compare two variances
##
## data: expend by stature
## F = 0.78445, num df = 12, denom df = 8, p-value = 0.6797
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1867876 2.7547991
## sample estimates:
## ratio of variances
## 0.784446
这里检验结果不显著,所以不能拒绝方差相等这个假设,但是置信区间非常宽。
5.5 两样本Wilcoxon检验
如果对正态分布假设有所怀疑,那么可能更需要使用一个非参数检验。两样本Wilcoxon检验用数据的秩(不考虑分组)代替数据本身,然后计算某一组中的秩和。这样便简化成了从1~ n1+n2中不重复地抽出n1个数字的问题。
wilcox.test(expend~stature)
## Warning in wilcox.test.default(x = c(7.53, 7.48, 8.08, 8.09, 10.15, 8.4, :
## cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: expend by stature
## W = 12, p-value = 0.002122
## alternative hypothesis: true location shift is not equal to 0
统计量 W 是第一组数据的秩和减去理论最小值(如果最小的n1个数都在第一组中,那么 W 就等于0)。
5.6 配对 t 检验
在同一个实验单位中有着两个度量值时使用配对检验。该检验主要通过做差来将问题简化为单样本检验。不过这种方法意味着我们要假设这个差值与不同水平的度量值是独立的。我们可以将每一对数构成的点与直线y=x画在同一幅图上,或是将每一对数的差与它们的均值画在同一幅图上(即Bland~Altman图),这是有效的图形检查方法。
#还是月经前后能量摄入
attach(intake)
intake
## pre post
## 1 5260 3910
## 2 5470 4220
## 3 5640 3885
## 4 6180 5160
## 5 6390 5645
## 6 6515 4680
## 7 6805 5265
## 8 7515 5975
## 9 7515 6790
## 10 8230 6900
## 11 8770 7335
#这里的关键是11位女性都被测量了2次,所以查看个人数据的差值是合理的
post-pre
## [1] -1350 -1250 -1755 -1020 -745 -1835 -1540 -1540 -725 -1330 -1435
#基本上都是负数
#这意味着月经后有更低的能量摄入
#t检验中paired=T表示希望进行一个配对检验
t.test(pre, post, paired = T)
##
## Paired t-test
##
## data: pre and post
## t = 11.941, df = 10, p-value = 3.059e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1074.072 1566.838
## sample estimates:
## mean of the differences
## 1320.455
t.test(pre, post)
##
## Welch Two Sample t-test
##
## data: pre and post
## t = 2.6242, df = 19.92, p-value = 0.01629
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 270.5633 2370.3458
## sample estimates:
## mean of x mean of y
## 6753.636 5433.182
#后者的t明显变小了,虽然在5%的水平下显著
#置信区间几乎是配对检验中的4倍
#从这两个可以看出,如果没有对比同一个人测量的前后信息,那么准确性就会降低
#当然也可以得出这样的结论
#在同一个人上进行两次测量,比对两组分别处于月经前、后的女性进行测量,效率更高
5.7 配对Wilcoxon检验
#与t.test函数相同
wilcox.test(pre, post, paired = T)
## Warning in wilcox.test.default(pre, post, paired = T): cannot compute exact
## p-value with ties
##
## Wilcoxon signed rank test with continuity correction
##
## data: pre and post
## V = 66, p-value = 0.00384
## alternative hypothesis: true location shift is not equal to 0
这个结果与 t 检验没有多大差别。p 值没有那么极端。
又一次,在由数据相等的情况下无法精确的计算 p 值,这里有两个相等的差值都是-1540。
对当前数据很容易计算精确的 p 值。这是11个正差值的概率,加上11个负差值的概率,即\(2\times(\frac{1}{2})^{11}=\frac{1}{1024}=0.00098\),所以近似计算的p 值差不多是真实值的4倍。
练习题
- 数据集react(注意这是一个向量)的值看起来是正态分布的吗,它的均值在 t 检验下显著地不等于0吗
- 在数据集vitcap中使用 t 检验比较两组肺活量,并计算99%置信区间。这个比较结果可能会产生误导,为什么
- 用非参数方法对react和vitcap两个数据做分析。
- 使用图像法检查intake数据集的配对 t 检验的建设是否合理。
- 函数shapiro.test基于Q-Q图的线性性来检验正态性。对react数据进行这个检验,它对移除异常值有帮助吗。
- 如果忽略潜在的时间因素影响,能通过简单的方法(怎么做)送ashina里的交叉实验中分析药物效果。提示,考虑个体内部的差异,如果只出现了时间因素的影响,两组差异会有什么表现
- 在10个含有25个观测值的模拟正态分布数据上分别做单样本 t 检验。用不同的分布再重复这个实验;尝试自由度为2的t 分布,以及指数分布(在后一种情况下对均值为1的分布做检验)。你能自动化这个实验从而重复更多次吗。