数据准备
df <- read.table(file = "D:/Documents/R wd/df.csv", header = T, sep = ",") # 数据导入。
df # 查看数据。
df$nitrogen <- as.factor(df$nitrogen) # 转化因子。
df$variety <- as.factor(df$variety) # 转化因子。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
第12章 重抽样与自助法
许多实际情况中统计假设（假定观测数据抽样自正态分布或者其他性质较好的理论分布）并不一定满足，比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况，这时基于随机化和重抽样的统计方法就可派上用场。
12.1 置换检验
置换检验的定义
置换检验(Permutation test)，也称随机化检验或重随机化检验，是Fisher于20世纪30年代提出的一种基于大量计算（computationally intensive），利用样本数据的全（或随机）排列，进行统计推断的方法，因其对总体分布自由，应用较为广泛，特别适用于总体分布未知的小样本资料，以及某些难以用常规方法分析资料的假设检验问题。
置换检验的原理
1、提出原假设，比如XX处理后结果没有变化
2、计算统计量，如两组的均值之差，记作t0
3、将所有样本放在一起，然后随机排序进行分组，再计算其统计量t1
4、重复第3步骤，直至所有排序可能性都齐全（比如有A组有n样本，B组有m样本，则总重复次数相当于从n+m中随机抽取n个的次数），得到一系列的统计量（t1-tn）
5、最后将这些统计量按照从小到大排序，构成抽样分布，再看t0是否落在分布的置信区间内（如95%置信区间），这时候可计算一个P值（如果抽样总体1000次统计量中大于t0的有10个，则估计的P值为10/1000=0.01），落在置信区间外则拒绝原假设
6、如果第3步骤是将所有可能性都计算了的话，则是精确检验；如果只取了计算了部分组合，则是近似结果，这时一般用蒙特卡罗模拟（Monte Carlo simulation）的方法进行置换检验
7、置换检验和参数检验都计算了统计量，但是前者是跟置换观测数据后获得的经验分布进行比较，后者则是跟理论分布进行比较。
请牢记：置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样（当做近似检验时）。因此，每次检验的结果都有所不同。
12.2 用coin包做置换检验
coin包提供了一个进行置换检验的一般性框架。通过该包，你可以回答如下问题。
响应值与组的分配独立吗？
两个数值变量独立吗？
两个类别型变量独立吗？
image.png
表12-2列出来的每个函数都是如下形式：
function_name(formula, data, distribution=)
其中：
formula描述的是要检验变量间的关系。示例可参见表12-2；
data是一个数据框；
distribution指定经验分布在零假设条件下的形式，可能值有exact，asymptotic和
approximate。
若distribution = "exact"，那么在零假设条件下，分布的计算是精确的（即依据所有可能的排列组合）。当然，也可以根据它的渐进分布（distribution = "asymptotic"）或蒙特卡洛重抽样（distribution = "approxiamate(B = #)"）来做近似计算，其中#指所需重复的次数。
distribution = "exact"当前仅可用于两样本问题。
12.2.1 独立两样本和K样本检验
library(coin) # 调用coin包。
df9 <- data.frame(score=c(40,57,45,55,58,57,64,55,62,65),treatment=factor(c(rep("A",5),rep("B",5)))) # 构建新数据集df9。
df9 # 显示数据。
## score treatment
## 1 40 A
## 2 57 A
## 3 45 A
## 4 55 A
## 5 58 A
## 6 57 B
## 7 64 B
## 8 55 B
## 9 62 B
## 10 65 B
t.test(score~treatment,data = df9, var.equal = TRUE) # t检验。
##
## Two Sample t-test
##
## data: score by treatment
## t = -2.345, df = 8, p-value = 0.04705
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
## -19.0405455 -0.1594545
## sample estimates:
## mean in group A mean in group B
## 51.0 60.6
oneway_test(score~treatment,data = df9, distribution = "exact") # 置换检验。
##
## Exact Two-Sample Fisher-Pitman Permutation Test
##
## data: score by treatment (A, B)
## Z = -1.9147, p-value = 0.07143
## alternative hypothesis: true mu is not equal to 0
传统t检验表明存在显著性差异（p < 0.05），而精确检验却表明差异并不显著（p > 0.072）。
第7章我用自己的数据进行了t检验，对比一下传统t检验和置换检验，结果如下：
t.test(v1 ~ nitrogen, data = df) # 传统t检验。
##
## Welch Two Sample t-test
##
## data: v1 by nitrogen
## t = -3.2071, df = 21.312, p-value = 0.004178
## alternative hypothesis: true difference in means between group N1 and group N2 is not equal to 0
## 95 percent confidence interval:
## -0.18538442 -0.03961558
## sample estimates:
## mean in group N1 mean in group N2
## 1.2000 1.3125
library(coin) # 调用coin包。
oneway_test(v1 ~ nitrogen, data = df, distribution = "exact") # 置换检验。
##
## Exact Two-Sample Fisher-Pitman Permutation Test
##
## data: v1 by nitrogen (N1, N2)
## Z = -2.7069, p-value = 0.004688
## alternative hypothesis: true mu is not equal to 0
两种检验方式下结果都是显著的
Wilcoxon-Mann-Whitney U检验
df10 <- data.frame(score1=c(10,27,25,155,158,166,187,174,155,22,25,22),trt1=factor(c(rep("A", 6),rep("B", 6)))) # 构建新数据集df10。
df10 # 显示数据。
## score1 trt1
## 1 10 A
## 2 27 A
## 3 25 A
## 4 155 A
## 5 158 A
## 6 166 A
## 7 187 B
## 8 174 B
## 9 155 B
## 10 22 B
## 11 25 B
## 12 22 B
with(df10, shapiro.test(score1[trt1 == "A"])) # 分类变量A对应的因变量score1数据正态性检验。
##
## Shapiro-Wilk normality test
##
## data: score1[trt1 == "A"]
## W = 0.75915, p-value = 0.02445
with(df10, shapiro.test(score1[trt1 == "B"])) # 分类变量B对应的因变量score1数据正态性检验
##
## Shapiro-Wilk normality test
##
## data: score1[trt1 == "B"]
## W = 0.75902, p-value = 0.02437
wilcox.test(score1 ~ trt1, data = df10, paired = TRUE) # 传统检验。
##
## Wilcoxon signed rank test with continuity correction
##
## data: score1 by trt1
## V = 9, p-value = 0.8335
## alternative hypothesis: true location shift is not equal to 0
library(coin) # 调用coin包。
wilcox_test(score1 ~ trt1, data = df10, distribution = "exact") # 置换检验。
##
## Exact Wilcoxon-Mann-Whitney Test
##
## data: score1 by trt1 (A, B)
## Z = -0.16097, p-value = 0.9134
## alternative hypothesis: true mu is not equal to 0
coin包规定所有的类别型变量都必须以因子形式编码。
wilcox.test()默认计算的也是精确分布。
K样本检验的置换检验
library(multcomp) # 调用multcomp包。
set.seed(1234) # 设置种子1234。
oneway_test(response~trt, data=cholesterol, distribution=approximate(nresample=9999)) # K样本检验。
##
## Approximative K-Sample Fisher-Pitman Permutation Test
##
## data: response by
## trt (1time, 2times, 4times, drugD, drugE)
## chi-squared = 36.381, p-value < 1e-04
12.2.2 列联表中的独立性
通过chisq_test()或cmh_test()函数，我们可用置换检验判断两类别型变量的独立性。 当数据可根据第三个类别型变量进行分层时，需要使用后一个函数。若变量都是有序型，可使用lbl_test()函数来检验是否存在线性趋势。
卡方独立性检验
library(vcd) # 调用vcd包。
mytable8 <- xtabs(~ Treatment + Improved, data = Arthritis) # 构建treatment和improved的列联表。
chisq.test(mytable8) # 卡方独立性检验。
##
## Pearson's Chi-squared test
##
## data: mytable8
## X-squared = 13.055, df = 2, p-value = 0.001463
卡方独立性检验的置换检验
library(coin) # 调用coin包。
library(vcd) # 调用vcd包。
Arthritis <- transform(Arthritis,
Improved = as.factor(as.numeric(Improved))) # 设置数据集中improved为因子。
set.seed(1234) # 设置种子。
chisq_test(Treatment~Improved, data=Arthritis,
distribution=approximate(nresample=9999)) # 卡方独立性检验的置换检验。
##
## Approximative Pearson Chi-Squared Test
##
## data: Treatment by Improved (1, 2, 3)
## chi-squared = 13.055, p-value = 0.0018
你可能会有疑问，为什么需要把变量Improved从一个有序因子变成一个分类因子？（好问题！）这是因为，如果你用有序因子，coin()将会生成一个线性与线性趋势检验，而不是卡方检验。
结果解读：两种检验下p值都是小于0.05，说明Treatment和Improved之间相互不独立
自己数据的演示
mytable3 <- xtabs(~ variety + nitrogen, data = df) # 构建列联表。
chisq.test(mytable3) # 传统因子间独立性检验。
##
## Pearson's Chi-squared test
##
## data: mytable3
## X-squared = 0, df = 1, p-value = 1
library(coin) # 调用coin包。
chisq_test(nitrogen ~ variety, data=df,distribution=approximate(nresample=9999)) # 独立性检验的置换检验。
##
## Approximative Pearson Chi-Squared Test
##
## data: nitrogen by variety (a, b)
## chi-squared = 0, p-value = 1
结果解读：p值均为1，表明nitrogen和variety相互独立。
12.2.3 数值变量间的独立性
spearman_test()函数提供了两数值变量的独立性置换检验。
states <- as.data.frame(state.x77) # 构建数据集states。
set.seed(1234) # 设置种子。
spearman_test(Illiteracy ~ Murder, data=states,
distribution=approximate(nresample=9999)) # 相关性置换检验。
##
## Approximative Spearman Correlation Test
##
## data: Illiteracy by Murder
## Z = 4.7065, p-value < 1e-04
## alternative hypothesis: true rho is not equal to 0
spearman_test(v1 ~ v2, data=df, distribution=approximate(nresample=9999)) # p值小于0.05，表明两个变量间独立性假设不满足。
##
## Approximative Spearman Correlation Test
##
## data: v1 by v2
## Z = 2.9455, p-value = 0.0016
## alternative hypothesis: true rho is not equal to 0
12.2.4 两样本和K样本相关性检验
当处于不同组的观测已经被分配得当，或者使用了重复测量时，样本相关检验便可派上用场。
对于两配对组的置换检验，可使用wilcoxsign_test()函数；多于两组时，使用friedman_test()函数。
library(coin) # 调用coin包。
library(MASS) # 调用MASS包。
wilcoxsign_test(U1 ~ U2, data=UScrime, distribution="exact") # 样本相关性置换检验。
##
## Exact Wilcoxon-Pratt Signed-Rank Test
##
## data: y by x (pos, neg)
## stratified by block
## Z = 5.9691, p-value = 1.421e-14
## alternative hypothesis: true mu is not equal to 0
自己数据演示
df10 <- data.frame(N=rep(c("N1","N2"),c(5,5)),y=c(1,2,3,4,5,11,12,13,14,15)) # 构建一个数据集df10，y1，y2分别是两个氮水平下的产量，各测定了5次。
t.test(y~N,data=df10) # 传统检验。结果表明N1和N2条件下产量存在差异。
##
## Welch Two Sample t-test
##
## data: y by N
## t = -10, df = 8, p-value = 8.488e-06
## alternative hypothesis: true difference in means between group N1 and group N2 is not equal to 0
## 95 percent confidence interval:
## -12.306004 -7.693996
## sample estimates:
## mean in group N1 mean in group N2
## 3 13
library(coin)
library(MASS)
wilcoxsign_test(df10$y[1:5] ~ df10$y[6:10], data=df10, distribution="exact") # 精确检验，结果表明N1和N2条件下产量存在差异。
##
## Exact Wilcoxon-Pratt Signed-Rank Test
##
## data: y by x (pos, neg)
## stratified by block
## Z = -2.2361, p-value = 0.0625
## alternative hypothesis: true mu is not equal to 0
12.2.5 深入探究
12.3 lmPerm包做置换检验
lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版，能够进行置换检验，而非正态理论检验。
lmp()和aovp()函数的参数与lm()和aov()函数类似，只额外添加了perm =参数。
perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样，直至估计的标准差在估计的p值0.1之下，判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意，若观测数大于10，perm = "Exact"将自动默认转为perm = "Prob"，因为精确检验只适用于小样本问题。
12.3.1 简单回归和多项式回归
简单线性回归的置换检验
R语言实战的例子：
library(lmPerm)
set.seed(1234)
fit19 <- lmp(weight ~ height, data=women, perm="Prob")
summary(fit19)
##
## Call:
## lmp(formula = weight ~ height, data = women, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7333 -1.1333 -0.3833 0.7417 3.1167
##
## Coefficients:
## Estimate Iter Pr(Prob)
## height 3.45 5000 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-Squared: 0.991, Adjusted R-squared: 0.9903
## F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
summary(lm(v1 ~ v3, data = df)) # 简单线性回归。
##
## Call:
## lm(formula = v1 ~ v3, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20656 -0.05668 0.02228 0.06561 0.19328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.253721 0.047939 26.153 <2e-16 ***
## v3 0.001667 0.028334 0.059 0.954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1041 on 22 degrees of freedom
## Multiple R-squared: 0.0001574, Adjusted R-squared: -0.04529
## F-statistic: 0.003463 on 1 and 22 DF, p-value: 0.9536
library(lmPerm) # 调用lmPerm包。
summary(lmp(v1 ~ v3, data = df, perm="Prob")) # 简单线性回归的置换检验。
##
## Call:
## lmp(formula = v1 ~ v3, data = df, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20656 -0.05668 0.02228 0.06561 0.19328
##
## Coefficients:
## Estimate Iter Pr(Prob)
## v3 0.001667 51 1
##
## Residual standard error: 0.1041 on 22 degrees of freedom
## Multiple R-Squared: 0.0001574, Adjusted R-squared: -0.04529
## F-statistic: 0.003463 on 1 and 22 DF, p-value: 0.9536
多项式回归的置换检验
R语言实战的例子：
library(lmPerm)
set.seed(1234)
fit20 <- lmp(weight ~ height + I(height^2), data=women, perm="Prob")
summary(fit20)
##
## Call:
## lmp(formula = weight ~ height + I(height^2), data = women, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.509405 -0.296105 -0.009405 0.286151 0.597059
##
## Coefficients:
## Estimate Iter Pr(Prob)
## height -7.34832 5000 <2e-16 ***
## I(height^2) 0.08306 5000 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-Squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
自己数据集的例子：
summary(lm(v1 ~ v3 + I(v3^2), data = df)) # 多项式回归。
##
## Call:
## lm(formula = v1 ~ v3 + I(v3^2), data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21079 -0.04893 0.01915 0.06477 0.18953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.240669 0.074300 16.698 1.34e-13 ***
## v3 0.026036 0.108218 0.241 0.812
## I(v3^2) -0.008352 0.035736 -0.234 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1064 on 21 degrees of freedom
## Multiple R-squared: 0.002751, Adjusted R-squared: -0.09223
## F-statistic: 0.02897 on 2 and 21 DF, p-value: 0.9715
library(lmPerm) # 调用lmPerm包。
summary(lmp(v1 ~ v3 + I(v3^2), data = df, perm="Prob")) # 多项式回归的置换检验。
##
## Call:
## lmp(formula = v1 ~ v3 + I(v3^2), data = df, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21079 -0.04893 0.01915 0.06477 0.18953
##
## Coefficients:
## Estimate Iter Pr(Prob)
## v3 0.026036 51 1.000
## I(v3^2) -0.008352 51 0.863
##
## Residual standard error: 0.1064 on 21 degrees of freedom
## Multiple R-Squared: 0.002751, Adjusted R-squared: -0.09223
## F-statistic: 0.02897 on 2 and 21 DF, p-value: 0.9715
12.3.2 多元回归
R语言实战的例子：
set.seed(1234)
states <- as.data.frame(state.x77)
fit21 <- lmp(Murder ~ Population + Illiteracy+Income+Frost,data=states, perm="Prob")
summary(fit21)
##
## Call:
## lmp(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.79597 -1.64946 -0.08112 1.48150 7.62104
##
## Coefficients:
## Estimate Iter Pr(Prob)
## Population 2.237e-04 51 1.0000
## Illiteracy 4.143e+00 5000 0.0004 ***
## Income 6.442e-05 51 1.0000
## Frost 5.813e-04 51 0.8627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-Squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
自己数据集的例子：
summary(lm(v1 ~ v3 + v4 + v2, data = df)) # 多元回归。
##
## Call:
## lm(formula = v1 ~ v3 + v4 + v2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16975 -0.04108 0.00548 0.04266 0.13547
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.88370 0.08700 10.158 2.43e-09 ***
## v3 0.03800 0.02600 1.461 0.1594
## v4 0.07547 0.02872 2.628 0.0161 *
## v2 -0.02210 0.03979 -0.556 0.5847
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07426 on 20 degrees of freedom
## Multiple R-squared: 0.5373, Adjusted R-squared: 0.4679
## F-statistic: 7.743 on 3 and 20 DF, p-value: 0.001265
library(lmPerm) # 调用lmPerm包。
summary(lmp(v1 ~ v3 + v4 + v2, data = df, perm="Prob")) # 多元回归的置换检验。
##
## Call:
## lmp(formula = v1 ~ v3 + v4 + v2, data = df, perm = "Prob")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16975 -0.04108 0.00548 0.04266 0.13547
##
## Coefficients:
## Estimate Iter Pr(Prob)
## v3 0.03800 469 0.177
## v4 0.07547 5000 0.012 *
## v2 -0.02210 71 0.592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07426 on 20 degrees of freedom
## Multiple R-Squared: 0.5373, Adjusted R-squared: 0.4679
## F-statistic: 7.743 on 3 and 20 DF, p-value: 0.001265
当两种方法所得结果不一致时，你需要更加谨慎地审视数据，这很可能是因为违反了正态性假设或者存在离群点。
12.3.3 单因素方差分析和协方差分析
R语言实战的例子：
library(lmPerm)
library(multcomp)
set.seed(1234)
fit22 <- aovp(response ~ trt, data=cholesterol, perm="Prob")
anova(fit22)
## Analysis of Variance Table
##
## Response: response
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## trt 4 1351.37 337.84 5000 < 2.2e-16 ***
## Residuals 45 468.75 10.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
自己数据集的例子：
summary(aov(v1 ~ nitrogen, data = df)) # 单因素方差分析。
## Df Sum Sq Mean Sq F value Pr(>F)
## nitrogen 1 0.07594 0.07594 10.29 0.00406 **
## Residuals 22 0.16243 0.00738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(lmPerm) # 调用lmPerm包。
summary(aovp(v1 ~ nitrogen, data = df, perm="Prob")) # 单因素方差分析的置换检验。
## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## nitrogen 1 0.075937 0.075937 5000 0.0022 **
## Residuals 22 0.162425 0.007383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R语言实战的例子：
library(lmPerm)
set.seed(1234)
fit23 <- lmp(weight ~ gesttime + dose, data=litter, perm="Prob")
anova(fit23)
## Analysis of Variance Table
##
## Response: weight
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## gesttime 1 161.49 161.493 5000 0.0006 ***
## dose 3 137.12 45.708 5000 0.0392 *
## Residuals 69 1151.27 16.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
自己数据集的例子：
summary(aov(v1 ~ nitrogen + block, data = df)) # 协方差分析。
## Df Sum Sq Mean Sq F value Pr(>F)
## nitrogen 1 0.07594 0.07594 10.170 0.00441 **
## block 1 0.00562 0.00562 0.753 0.39523
## Residuals 21 0.15680 0.00747
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(lmPerm) # 调用lmPerm包。
summary(aovp(v1 ~ nitrogen + block, data = df, perm="Prob")) # 协方差分析的置换检验。
## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## nitrogen 1 0.075937 0.075937 5000 0.0018 **
## block 1 0.005625 0.005625 185 0.3514
## Residuals 21 0.156800 0.007467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
12.3.4 双因素方差分析
R语言实战的例子：
library(lmPerm)
set.seed(1234)
fit24 <- lmp(len ~ supp*dose, data=ToothGrowth, perm="Prob")
anova(fit24)
## Analysis of Variance Table
##
## Response: len
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## supp 1 205.35 205.35 5000 < 2e-16 ***
## dose 1 2224.30 2224.30 5000 < 2e-16 ***
## supp:dose 1 88.92 88.92 2032 0.04724 *
## Residuals 56 933.63 16.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
自己数据集的例子：
summary(aov(v1 ~ year*nitrogen, data = df)) # 双因素方差分析。
## Df Sum Sq Mean Sq F value Pr(>F)
## year 1 0.00510 0.00510 0.664 0.42486
## nitrogen 1 0.07594 0.07594 9.874 0.00513 **
## year:nitrogen 1 0.00350 0.00350 0.456 0.50740
## Residuals 20 0.15382 0.00769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(lmPerm) # 调用lmPerm包。
summary(aovp(v1 ~ year*nitrogen, data = df, perm="Prob")) # 双因素方差分析的置换检验。
## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## year 1 0.005104 0.005104 143 0.4126
## nitrogen 1 0.075937 0.075937 5000 0.0050 **
## year:nitrogen 1 0.003504 0.003504 51 0.7255
## Residuals 20 0.153817 0.007691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
值得注意的是，当将aovp()应用到方差分析设计中时，它默认使用唯一平方和法（SAS也称为类型III平方和）。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和（SAS是类型I平方和）。每种效应依据模型中先出现的效应做相应调整。对于平衡设计，两种方法结果相同，但是对于每个单元格观测数不同的不平衡设计，两种方法结果则不同。不平衡性越大，结果分歧越大。若在aovp()函数中设定seqs = TRUE，可以生成你想要的序贯平方和。
12.4 置换检验点评
你可能已经注意到，基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好，两种方法结果的一致性也验证了正态理论方法适用于上述示例。
当然，置换检验真正发挥功用的地方是处理非正态数据（如分布偏倚很大）、存在离群点、样本很小或无法做参数检验等情况。不过，如果初始样本对感兴趣的总体情况代表性很差，即使是置换检验也无法提高推断效果。
置换检验主要用于生成检验零假设的p值，它有助于回答“效应是否存在”这样的问题。不过，置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是，这正是自助法大显神通的地方。
12.5 自助法
所谓自助法，即从初始样本重复随机替换抽样，生成一个或一系列待检验统计量的经验分布。 无需假设一个特定的理论分布，便可生成统计量的置信区间，并能检验统计假设。
倘若你假设均值的样本分布不是正态分布，该怎么办呢？可使用自助法。
(1)从样本中随机选择10个观测，抽样后再放回。有些观测可能会被选择多次，有些可能一直都不会被选中。
(2)计算并记录样本均值。
(3)重复1和2一千次。
(4)将1000个样本均值从小到大排序。
(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数，它们就限定了95%的置信区间。
12.6 boot包中的自助法
boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。
一般来说，自助法有三个主要步骤。
(1)写一个能返回待研究统计量值的函数。如果只有单个统计量（如中位数），函数应该返回一个数值；如果有一列统计量（如一列回归系数），函数应该返回一个向量。
(2)为生成R中自助法所需的有效统计量重复数，使用boot()函数对上面所写的函数进行处理。
(3)使用boot.ci()函数获取第(2)步生成的统计量的置信区间。
主要的自助法函数是boot()，它的格式为：
bootobject <- boot(data=, statistic=, R=, ...)
参数见下表：
image.png
boot()函数调用统计量函数R次，每次都从整数1:nrow(data)中生成一列有放回的随机指标，这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算，结果存储在bootobject中。
image.png
你可以用bootobject t0和bootobject t来获取这些元素。
一旦生成了自助样本，可通过print()和plot()来检查结果。如果结果看起来还算合理， 使用boot.ci()函数获取统计量的置信区间。格式如下:
boot.ci(bootobject, conf=, type= )
type参数设定了获取置信区间的方法。perc方法（分位数）展示的是样本均值，bca将根据偏差对区间做简单调整。
12.6.1 对单个统计量使用自助法
回归的R平方值
rsq <- function(formula, data, indices) {
d <- data[indices,]
fit25 <- lm(formula, data=d)
return(summary(fit25)$r.square)
}
1000次自助抽样
library(boot)
set.seed(1234)
results1 <- boot(data=mtcars, statistic=rsq,
R=1000, formula=mpg~wt+disp)
输出结果
print(results1)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~
## wt + disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7809306 0.01379126 0.05113904
结果可视化
plot(results1)
image.png
95%的置信区间获取
boot.ci(results1, type=c("perc", "bca"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results1, type = c("perc", "bca"))
##
## Intervals :
## Level Percentile BCa
## 95% ( 0.6753, 0.8835 ) ( 0.6344, 0.8561 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
12.6.2 多个统计量的自助法
回归系数向量函数
bs <- function(formula, data, indices) {
d <- data[indices,]
fit26 <- lm(formula, data=d)
return(coef(fit26))
}
自助抽样1000次
library(boot)
set.seed(1234)
results2 <- boot(data=mtcars, statistic=bs,
R=1000, formula=mpg~wt+disp)
print(results2)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mtcars, statistic = bs, R = 1000, formula = mpg ~
## wt + disp)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 34.96055404 4.715497e-02 2.546106756
## t2* -3.35082533 -4.908125e-02 1.154800744
## t3* -0.01772474 6.230927e-05 0.008518022
获得车重和发动机排量95%的置信区间
boot.ci(results2, type="bca", index=2)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results2, type = "bca", index = 2)
##
## Intervals :
## Level BCa
## 95% (-5.477, -0.937 )
## Calculations and Intervals on Original Scale
boot.ci(results2, type="bca", index=3)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = results2, type = "bca", index = 3)
##
## Intervals :
## Level BCa
## 95% (-0.0334, -0.0011 )
## Calculations and Intervals on Original Scale
12.7 小结
置换检验和自助法并不是万能的，它们无法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳，或者样本量过小而无法准确地反映总体情况，这些方法也是爱莫能助。
参考资料：
- 《R语言实战》（中文版），人民邮电出版社，2013.
- 置换检验（Permutation Test），https://www.bioinfo-scrounger.com/archives/564/
- 统计上的“置换检验”有什么用？，https://www.zhihu.com/question/26629875/answer/392615941
