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)