第12章 重抽样与自助法

R语言的再复习之路
 
 

置换检验

置换检验步骤：
（1）与参数方法类似，计算观测数据的t统计量，称为t0；
（2）将两组数据放在一个组中；
（3）随机分配一半到A处理中，分配一半到B处理中；
（4）计算并记录新观测的t统计量；
（5）对每一种可能的随机分配重复步骤（3）~（4）；
（6）将所有情况下的t统计量按升序排列，这便是基于样本数据的经验分布；
（7）如果t0落在经验分布中间95%部分的外面，则在0.05的显著性水平下，拒绝两个处理组的总体均值相等的零假设。
 
注：置换方法与参数方法都计算了相同的t统计量，区别在于参数方法是将统计量与理论分布进行比较，而置换方法则是将t统计量与置换观测数据后获得的经验分布进行比较。
 

1. coin包

检验	coin函数
两样本和K样本置换检验	oneway_test(y ~ A)
含一个分层（区组）因子的两样本和K样本置换检验	oneway_test(y ~ A | C)
Wilcoxon-Mann-Whitney秩和检验	wilcox_test(y ~ A)
Kruskal-Wallis检验	kruskal_test(y ~ A)
Pearson卡方检验	chisq_test(A ~ B)
Cochran-Mantel-Haenszel检验	cmh_test(A ~ B | C)
线性相关检验	lbl_test(D ~ E)
Spearman检验	spearman_test(y ~ x)
Friedman检验	friedman_test(y ~ A | C)
Wilcoxon秩和检验	wilcoxonsign_test(y1 ~ y2)

格式：function_name(formula, data, distribution = )，其中distribution指定经验分布在零假设条件下的形式，可能值有exact、asymptotic和approxmate。

1.1 独立两样本和K样本检验

# 创建数据
library(coin)
score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)
treatment <- factor(c(rep('A', 5), rep('B', 5)))
mydata <- data.frame(treatment, score)

# 参数方法
t.test(score ~ treatment, data = mydata, var.equal = TRUE)

# 置换方法
oneway_test(score ~ treatment, data = mydata, distribution = 'exact')

# K样本检验
library(multcomp)
set.seed(1234)
oneway_test(response ~ trt, data = cholesterol, distribution = approximate(B = 9999))

1.2 列联表中的独立性

library(coin)
library(vcd)
Arthritis <- transform(Arthritis, Improved = as.factor(as.numeric(Improved)))
set.seed(1234)
chisq_test(Treatment ~ Improved, data = Arthritis, distribution = approximate(B = 9999))

注：这里把变量Improved从一个有序因子变为一个分类因子。因为如果用有序因子，coin()函数将会生成一个线性与线性趋势检验，而不是卡方检验。

1.3 数值变量间的独立性

states <- as.data.frame(state.x77)
set.seed(1234)
spearman_test(Illiteracy ~ Murder, data = states, distribution = approximate(B = 9999))

1.4 两样本和K样本相关性检验

对于两配对组的置换检验，可使用wilcoxsign_test()函数；多于两组时，使用friedman_test()函数。

library(coin)
library(MASS)
wilcoxsign_test(U1 ~ U2, data = UScrime, distribution = 'exact')

 

2. lmPerm包

 
 

自助法

自助法步骤：
（1）从样本中随机选择n个观测，抽样后再放回。有些观测可能会被选择多次，有些可能一直都不会被选中；
（2）计算并记录样本均值；
（3）重复（1）~（2）步骤1000次；
（4）将1000个样本均值从小到大排序；
（5）找出样本均值2.5%和97.5%的分位点，此时即初始位置和最末位置的第25个数，他们就限定了95%的置信区间。
 
格式：bootobject <- boot(data = , statistic = , R = , ...)
 

1. 对单个统计量使用自助法

# 获取R平方的函数
rsq <- function(formula, data, indices){
    d <- data[indices, ]
    fit <- lm(formula, data = d)
    return(summary(fit)$r.square)
    }
    
# 自助法抽样
library(boot)
set.seed(1234)
results <- boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~ wt + disp)

# 输出结果
print(results)
plot(results)
boot.ci(results, type = c('perc', 'bca'))

2. 多个统计量的自助法

# 获取模型参数的函数
bs <- function(formula, data, indices){
    d <- data[indices, ]
    fit <- lm(formula, data = d)
    return(coef(fit))

# 自助法抽样
library(boot)
set.seed(1234)
results <- boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~ wt + disp)

# 输出结果
print(results)
plot(results, index = 2)
boot.ci(results, type = 'bca', index = 2)