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)