r语言参数检验结果分析 R语言学习笔记九

导语：上一期跟大家介绍了R语言学习笔记（八） -极大似然估计，本期跟大家一起学习假设检验的基本原理和计算方法。

还是以抛硬币的例子来说明，假如我们想知道一枚硬币是否均匀，抛硬币100次后统计正反面向上的次数，得到59次正面和41次反面，明显偏离1:1，那么我们是否应该判断该硬币不均匀呢？

现在我们用R语言来实现这一推理过程，假设硬币是均匀的，那么正反面向上的概率各0.5，画出正面向上次数的概率分布：

> numFlips = 100 > numHeads = 59 > k = 0:numFlips > numHeads = sum(coinFlips == "H") > binomDensity = tibble(k = k, p = dbinom(k, size = numFlips, prob = 0.5)) > library("ggplot2") > ggplot(binomDensity) geom_bar(aes(x = k, y = p), stat = "identity") geom_vline(xintercept = numHeads, col = "blue"

从图中看到概率最大的k是50，也就是最有可能出现50次正面向上，这很容易理解。现在问题来了，我们的样本数据是59，出现50周围的数字的概率也比较大，我们到底应该判断硬币是否均匀呢？

抛100次硬币出现正面向上的次数k的范围是0-100，这里我们按概率将其分成两个区间，接受域和拒绝域。若出现的样本落在拒绝域中，则拒绝原假设，即认为硬币不均匀。我们还是用代码说明:

> library("dplyr") > alpha = 0.05 > binomDensity = arrange(binomDensity, p) %>% mutate(reject = (cumsum(p) <= alpha)) > ggplot(binomDensity) geom_bar(aes(x = k, y = p, col = reject), stat = "identity") scale_colour_manual( values = c(`TRUE` = "red", `FALSE` = "darkgrey")) geom_vline(xintercept = numHeads, col = "blue") theme(legend.position = "none")

解释一下上述代码。首先我们设定了拒绝域的比例α = 0.05（即显著性水平），然后用dplyr包中的arrange函数对p值从小到大排序。根据累积概率是否大于α将数据分为接受域和拒绝域。下图中红色部分即为拒绝域，我们的样本落在接受域，因此我们应该接受原假设，认为硬币就是均匀的。

事实上R有专门的函数完成上述过程，只需要一个简单的命令：

> binom.test(x = numHeads, n = numFlips, p = 0.5) Exact binomial test data: numHeads and numFlips number of successes = 59, number of trials = 100, p-value = 0.08863 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.4871442 0.6873800 sample estimates: probability of success 0.59

上述过程其实就是一个假设检验，p-value = 0.08863 > 0.05，因此应该接受原假设，即伯努利试验的成功概率等于0.5。

p-value的概念其实是比较难理解的，用通俗的话来说，就是在原假设（H0）正确时，出现现状或更极端的情况的概率。比如上面的例子中，p-value反映的是样本出现正面向上≥59或≤41的概率。我们需要记住的是，当p-value小于显著性水平α时，就拒绝原假设。

（1）根据实际问题的要求，提出原假设H0和备择假设H1；

（2）给定显著性水平α和样本容量n；

（3）确定检验统计量以及拒绝域的形式；

（4）按P{当H0为真拒绝H0} ≤ α求出拒绝域；

（5）根据样本观察值做出决策，接受或拒接H0。

我们再回顾一下假设检验的两类错误：

我们无法同时减少这两类错误，因为二者其实是此消彼长的关系，如下图所示，两个峰分别表示统计量在原假设和备择假设下的分布。假设决策边界是黑线，如果统计量位于黑线右侧，则拒绝原假设。假阳性（犯第Ⅰ类错误）的概率就是暗红色区域，假阴性（犯第Ⅱ类错误）的概率是深蓝色区域。当黑线右移时，假阳性减小但假阴性增大；黑线左移时则相反。因此更多的时候我们需要在两类错误中做一个权衡。

我们在实验中最常见的是需要对两组数据进行比较，比如测试某种药物对人体是否有效，某种肥料是否能使植物增产等等。这种两组数据的比较通常都需要进行t检验，下面还是用一个例子说明。PlantGrowth描述了植物在对照（ctrl）和两种处理（trt1和trt2）条件下的干重，我们先画出图像。

> library("ggbeeswarm") > data("PlantGrowth") > ggplot(PlantGrowth, aes(y = weight, x = group, col = group)) geom_beeswarm() theme(legend.position = "none")

现在我们想比较两种处理对植物干重的影响，做t检验如下：

> tt = with(PlantGrowth, t.test(weight[group =="ctrl"], weight[group =="trt2"], var.equal = TRUE)) > tt Two Sample t-test data:weight[group == "ctrl"] and weight[group == "trt2"] t = -2.134, df = 18, p-value = 0.04685 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.980338117 -0.007661883 sample estimates: mean of x mean of y 5.032 5.526

结果给最重要的信息有两个，一个是p-value = 0.04685 < 0.05，另外一个是95%置信区间 = (-0.980 ~ -0.008)不包含0。二者给出的结论是一致的，拒绝原假设，即认为两种处理对植物干重的影响是不一样的。

上述t检验的过程基于两个假设，首先是两个样本都来自正态总体，其次是样本方差相等，这其实是一种理想化的假设。我们可以用置换检验的方法对上述过程进行验证，方法是对样本数据进行重复抽样，做大量的t检验，得出一个统计量的分布。类似于蒙特卡洛方法，置换检验通常也需要很大的计算量。我们用R语言进行上述数据的置换检验：

> abs_t_null = with( dplyr::filter(PlantGrowth, group %in% c("ctrl", "trt2")), #用filter函数筛选 replicate(10000, abs(t.test(weight ~ sample(group))$statistic))) #用sample函数对样本标签随机置换 > ggplot(tibble(`|t|` = abs_t_null), aes(x = `|t|`)) geom_histogram(binwidth = 0.1, boundary = 0) geom_vline(xintercept = abs(tt$statistic), col = "red")

我们可以计算出样本统计量处在上述分布的哪个区间，也就是上图中红色线条右侧的面积占比。

> mean(abs(tt$statistic) <= abs_t_null) [1] 0.0489

这个数值跟我们上面t检验的p-value（0.04685）非常近似，因此可以认为上面的t测验的结果是合理的。