作者：姚某某
本节主要总结「数据分析」的「Resampling」重抽样思想，并通过 R 语言实现。
有一种东西叫作「传统」，它在很多时候很有用，但会让你思维固化，在新的环境下让你出错。
在总结回归分析和方差分析的时候 ④Ｒ语言之数据分析「初章」，我总是会在模型的建立之前提到「统计假设」，在模型建立之后进行「假设检验」，原因想必大家都能理解，就是因为这些「统计假设」是我们模型建立思想的基础，是支撑我们模型正确性的「必要条件」。但是，不可否认的是，这些「必要条件」最终会成为我们「数据分析」的局限，让我们对「不满足条件的数据集」束手无策。
本节，我们就来解决这个「必要条件」中的其中一条假设，从特例到普遍。
统计假设中有一条，叫做「假定观测数据抽样自正态分布或者是是其他性质较好的分布」，那么当数据集抽样自「未知分布、混合分布」、样本容量过小或存在离群点时，传统的统计方法所得到的模型可能就会不那么准确，原因之前已经讲过，这个时候「Resampling」 的思想就出现了。它抛弃了分布的理论，而是完全基于同一个样本，在这个样本中多次重复抽样，然后将所有抽样的结果作为总体，将原样本放到其中去检验，判定其显著性。因为需要多次重复抽样，所有它被称为重抽样「Resampling」。
1. Resampling 的分类
1.1. 置换检验（ permutation test ）
这个方法是传统统计方法的创建者 R. A. Fisher 建立的，但是由于这个方法的计算量过大、且计算机技术也未成熟，他最后放弃了这个方法。但是，数十年后的今天，计算机技术的高速发展，这个方法终于能够实现并发挥其价值。
为了更清楚的说明置换检验的思想，我举一个「两总体均值之差」的推断问题：
现在有两套学习方案 A 和 B，在 10 个受试者中随机抽取 5 个按照方案 A 学习，另外 5 个按照方案 B 学习，在学习完毕后对 10 个受试者进行测试，得到分数如下：
方案 A 方案 B
91 89
88 78
76 93
79 81
82 77
原假设H_{0} ：两种方案的总体均值相等 ；备择假设H_{a}：两种方案的总体均值不等
这种问题，用参数方法来解决你应该是很熟悉的，现在我来解释置换检验的思想是怎样的。置换检验的思想下，由于我们原假设两个方案的总体均值是相等的，那么如果两种方案真的是等价的话，则 10 个受试者的分数的分配应该是任意的。这么讲你可能还是不太理解，其实在原假设成立的前提下，我们把两种方案的总体视为等价的，那么我们同样也可将这两个方案的总体视为同一个总体，这样的话，10 个受试者者的分数在可以在任意位置，而不用在乎这个位置是属于方案 A 还是方案 B。在这样的前提下，我们利用重抽样得到一个「源于已有样本的抽样分布」，再利用 t 统计量在这个「源于已有样本的抽样分布」上判别初始观测数据的显著性，具体步骤如下：
1. 计算初始观测数据的 t 统计量，记为 t0
2. 将初始观测数据的 10 个得分，作为一个「重抽样的总体」，从中抽取 5 个放到 A 组，另五个放到 B 组，并计算此种情况下的 t 统计量
3. 重复步骤 2，直到穷尽所有可能。其实就是一个组合问题，一共有252种可能
4. 将 252 种可能的 t 统计量，按照从小到大排序，得到了「源于已有样本的抽样分布」
5. 根据t0在此分布中的位置，和我们确定的显著性水平，判别原假设我们可否拒绝
1.2. 交叉验证（ Cross validation ）
在交叉验证的思想中，一个样本被随机的分成两个或 K 个子样本，将其中一个作为验证集，其他 K-1 个子样本组合成训练集，在训练集上获得回归方程，而后在验证集上做出预测，这为一重验证，而后所有子样本轮流充当验证集，如此往复，直至遍历所有可能。这样我们就可以得到 K 个“预测误差”，求其平均值即为所谓的“ CV 值”，所以常说的 CV 值实际上是预测误差期望Err的一个估计值。
交叉验证的思想比较容易理解，也是一种十分符合我们直觉的验证方法，这里就不举例说明了。但值得注意的一点是，K 的取值会影响我们预测的准确性，具体是怎么影响的涉及到“偏误”、“方差”以及“计算效率”的权衡，主要是“偏误”与“计算效率”的权衡，这里就不深究，感兴趣的同学可以看看 「模型选择的一些基本思想和方法」这篇文章。不过，我们可以粗略的设置 K 的取值：
- 大样本，5 重交叉验证，对预测误差的估计已经足够准确，这里偏向于增大计算效率
- 小样本，10重交叉验证，为了得到相对准确的估计，这里考虑损失计算效率
1.3. 刀切法（ Jackknife ）
刀切法的思想，十分适用于分布的分散过宽或者出现异常值的数据，通俗来讲就是「既然我们所拥有的数据都只是样本，都是抽样得到的，那么我们在做推断和估计的时候，扔掉几个观测，看看效果如何。」在刀切法中我们反复的行为就是从样本中选取一个观测，把它抛弃，而后利用剩余的数据做出估计，如此往复，直到每一个观测都被抛弃过一次。根据每次估计的统计量，取均值，而后将这个均值与原样本（没有观测没抛弃）上估计到到的统计量想比较，就可以得到「偏差校正刀切估计值」。我的表述可能你能听懂，但是还理解不深，下面我举例说明。
因为一个理科生的各科成绩之间可能存在某种相关性，现在我们想用用学生的物理成绩和化学成绩来预测他的数学成绩，现在有 50 个观测：
序号 物理 化学 数学
1 78 87 90
2 88 79 83
3 67 71 78
… … … …
50 89 88 92
1.4. Bootstrap
与刀切法相比，Bootstrap 法的重抽样思想更加透彻。刀切法的重抽样次数与观测数量有关，而 Bootstrap 法的抽样次数理论上是没有上限的，将原样本在计算机性能允许的基础上尽可能重复抽样更多次，得到了一个「virtual population 」，从其中我们可以得到统计量的「虚拟潜在分布」，并根据此分布估计置信区间。
举例说明：
现在有一个容量为 10 的样本，均值 20，标准差 2，希望根据样本估计总体均值的置信区间。
按照传统思想，我们可以假设样本分布为正态分布，而后根据正态分布计算总体均值的置信区间。但当样本分布不服从正态分布时，可以利用 Bootstrap 法：1. 有放回地从样本中抽取 10 个观测，并计算其均值（注意由于抽样是有放回的，重抽样得到的样本与原样本虽然容量相同，但是观测不一定相同，因为有的观测可能被多次选择，有的可能一直都不会被选中）2. 重复 1，次数为 20003. 将 2000 个均值升序排列4. 确定显著性水平，计算得到置信区间（例如：要求 95% 的置信区间，找出 2.5% 和 97.5% 的分位点，其中间部分就为 95% 的置信区间）
2. R 语言实现重抽样
2.1. 置换检验
格式：
> function_name( formula , data , distribution = )
# formula 是被检验变量之间的关系， data 是数据框，distribution 指定重抽样的方式
这里的重抽样方式，除了之前将的按照完全排列组合抽样之外（ distribution = "exact" ），为了节省计算机资源，可以选择根据渐进分布抽样（ distribution = "asymptotic" ），和蒙特卡洛重抽样（ distribution = "approxiamate(B=#)" ）# 为重复次数，一般要设定随机种子，以便将来重现。
2.1.1. 两样本和 K 样本的独立性检验（ coin 包 ）
两样本独立性的参数检验法（ 置换检验 ）
> oneway_test(score~treatment,data=mydata,distribution="exact")# 这里表达式 ~ 的两边分别为应变量和分组变量# 传统为 t.test() 函数
两样本独立性的非参数检验法（ 置换检验 ）
> wilcox_test(Prob ~ So ,data=UScrime, distribution="exact")# 这里表达式 ~ 的两边分别为应变量和分组变量# 传统为 wilcox.test() 函数
K 样本独立性的参数检验法（ 置换检验 ）
> oneway_test(response ~ trt ,data=cholesterol,distribution="approxiamate(B=9999)")# 这里表达式 ~ 的两边分别为应变量和分组变量# 传统为 单因素方差分析
2.1.2. 列联表中两分组变量之间的独立性检验（ coin 包 ）
> Arthritis <- transform(Arthritis, Improved=as.factor(as.numeric(Improved)))> set.seed(1234)> chisq_test(Treatment~Improved,data=Arthritis, distribution=approximate(B=9999))# 这里表达式 ~ 的两边分别为两个分组变量# 如果使用有序因子，则次函数为线性趋势检验，因此第一行代码将 Improved 转换成无序
2.1.3. 数值型变量之间的独立性检验（ coin 包 ）
> states <- as.data.frame(state.x77)> set.seed(1234)> spearman_test(Illiteracy ~ Murder , data=states, distribution=approximate(B=9999))# 这里表达式 ~ 的两边分别为两个数值型变量# coin 包中必须为数据框
2.1.4. 两样本和 K 样本的相关性检验（ coin包 ）
两样本
> wilcoxsign_test(U1~U2,data=UScrime,distribution="exact")
K 样本
> friedman_test(times ~ methods | block, data = rounding)
2.1.5. 简单回归和多项式回归（ lmPerm 包 ）
简单线性回归（置换检验）
> set.seed(1234)> fit1 <- lmp(weight~height,data = women, perm="Prob")> summary(fit)# lmp() 函数与 lm() 函数类似，只是多了 perm 参数，其可选值有 Exact、Prob 或SPR# Exact 根据所有可能的排列组合生成精确检验。Prob 从所有可能的排列中不断抽样，直至估计的标准差在估计的 p 值 0.1 之下，判停准则由可选的 Ca 参数控制。SPR 使用贯序概率比检验来判断何时停止抽样。# 观测数大于10，perm="Exact"将 自动默认转为 perm="Prob"，因为精确检验只适用于小样本问题。
多项式回归（置换检验）
> set.seed(1234)> fit2 <- lmp(weight~height+I(height^2),data=women,perm="Prob")> summary(fit)
2.1.6. 多元回归（置换检验）
> set.seed(1234)> states <- as.data.frame(state.x77)> fit <- lmp(Murder ~ Population + Illiteracy + Income + Frost, data=states, perm = "Prob")> summary(fit)
2.1.7. 单因素方差分析和协方差分析
单因素方差分析（置换检验）
> set.seed(1234)> fit <- aovp(response ~ trt, data = cholesterol, perm="Prob")> summary(fit)# aovp() 函数与 aov() 函数类似，只是多了 perm 参数，其可选值有 Exact、Prob 或 SPR
单因素协方差分析（置换检验）
> set.seed(1234)> fit <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob")> summary(fit)
2.1.8. 双因素方差分析
> set.seed(1234)> fit <- aovp(len~supp*dose, data = ToothGrowth, perm="Prob")> summary(fit)
2.2. 交叉验证
《 R 语言实战 》中的自编函数 shrinkage ( ) 可以实现 K 重交叉验证，默认为 10 重，可以通过参数 k 修改。
> states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) > fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states) > shrinkage(fit)
2.3. 刀切法
《 R 语言实践 》中并未介绍刀切法的实现，可以参考以下代码
# 抽样得到样本> set.seed(1234)> x <- rnorm(100,0,5)> # 运用刀切法重抽样> jack <- function(x){+ jackknife <- 0+ for(i in 1:length(x)){+ jackknife[i]=length(x)*var(x)-(length(x)-1)/length(x)*sum(var(x[-i]))+ }+ return(jackknife)+ }> # 计算刀切法得到的抽样方差> mean(jack(x))/length(x)[1] 24.97106> # 计算抽样方差> var(x)[1] 25.22075
2.4. Bootstrap
格式：
# Bootstrap 重抽样，并返回统计量> bootobject <- boot(data=, statistic=, R=, ...)# data 为原样本，也是我们重抽样的数据来源，可以是向量、矩阵或者数据框# statistic 为生成 k 个统计量的函数，这个函数的参数中必须包括 indices，indices 参数用于 boot() 对原样本重抽样# R 为重抽样的次数# 获取置信区间> boot.ci(bootobject, conf=, type= )# bootobject 为存储重抽样结果的对象# conf 为置信区间的置信度# type 为置信区间的类型
2.4.1. 对单个统计量使用自助法
> #获取 R 平方值的自编函数> rsq <- function(formula, data, indices){+ d <- data[indices,]+ fit <- lm(formula, data=d)+ return(summary(fit)$r.square)+ }> #自助抽样> set.seed(1234)> results <- boot(data=mtcars, statistic = rsq,+ R=1000, formula=mpg~wt+disp)> #输出 boot 对象> print(results)> plot(results)> #计算置信区间> boot.ci(results, type=c("perc","bca"))# type 参数为 perc 时展示的是样本均值，为 bca 时将根据偏差对区间做简单调整
2.4.2. 多个统计量的自助法
> #返回回归系数向量的自编函数> bs <- function(formula, data, indices){+ d <- data[indices,]+ fit <- lm(formula,data=d)+ return(coef(fit))+ }> #自助抽样> set.seed(1234)> results <- boot(data=mtcars, statistic = bs,+ R=1000, formula=mpg~wt+disp)> #输出 boot 对象> print(results)> plot(results,index=2)> #计算置信区间> boot.ci(results,type = "bca", index=2)> boot.ci(results,type = "bca", index=3)# index 参数用于指明需要分析的重抽样所得样本的列