#重抽样与自助法
#置换检验
#又称为随机化检验,或重随机化检验
#假如你觉的假定数据成正态分布并不合适,或者担心离群点的影响,又或者感觉对于
#标准的参数方法来说数据集太小,那么置换检验便提供了一个非常不错的选择
#置换检验主要用于生成检验零假设的p值,它有助于回答"效应是否存在"这样的问题
#置换方法对于获取置信区间和估计测量精度是比较困难的
#对比传统基于正态理论的检验与置换检验
#结果非常接近
#置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大,存在离群点,样本很小
#或无法做参数检验等情况)
#置换检验的步骤
#1 与参数方法类似,计算观测数据的t统计量,称为t0;
#2 将10个得分放在一个组中
#3 随机分配五个得分到A处理,并分配五个得分到B处理
#4 计算并记录新观测的t统计量
#5 对每一种可能随机分配重复步骤3-4,此处有252中可能的分配组合
#6 将252个t统计量按升序排列,这便是基于(或以之为条件)样本数据的经验分布
#7 如果t0落在经验分布中间95%部分的外面,则在0.05的显著性水平下,拒绝两个处理组
#的总体均值相等的零假设
#结论:
#置换方法和参数方法都计算了相同的t统计量。
#但置换方法并不是将统计量与理论分布进行比较,而是将其与置换观测数据后获得的经验分布进行比较
#根据统计量值的极端性判断是否有足够的理由拒绝零假设。
#检验分布依据的是数据所有可能的排列组合,所以称为精确性检验
#coin包
#对于独立性问题提供了非常全面的置换检验的框架
#在coin包中,类别型变量和序数变量必须分别转化为因子和有序因子。另外数据要以数据框形式存储
#提供了一个置换检验的一般性框架,可以分析一组变量相对于其他任意变量,是否与第二组变量相互独立
#y和x是数值变量,A和B是分类因子,C是类别型区组变量,D和E是有序因子
#y1和y2是相匹配的数值变量
#distribution:指定经验分布在零假设条件下
#exact,分布计算是精确的,即依据所有可能的排列组合
#asymptotic,渐进分布
#approxiamate(B=#),蒙特卡洛重抽样,#表示所要重复的次数
install.packages(c("coin"))
library("survival")
library("coin")
#范例一
library(coin)
score <- c(40,57,45,55,58,57,64,55,62,65)
#a <- rep("A",5)
#class(a) #返回"character"型向量:"A" "A" "A" "A" "A"
##将字符长向量转化成因子长向量A A A A A B B B B B
treatment <- factor(c(rep("A",5),rep("B",5)))
mydata <- data.frame(treatment,score) #基于向量构造data.frame,向量=列向量
#传统t检验
t.test(score ~ treatment,data=mydata,var.equal=TRUE)
#p-value = 0.04705,p<0.05,落入拒绝域,拒绝相等的原假设,存在显著性差异
#置换检验
#oneway_test(y~A)
#两样本和K样本置换检验
oneway_test(score~treatment,data=mydata,distribution="exact")
#p-value = 0.07143,P>0.5,未落入拒绝域,接受相等的原假设,不存在显著性差异
#结论:
#由于只有10个观测,更倾向于相信置换检验的结果。
#范例二
library(MASS)
UScrime <- transform(UScrime,So=factor(So))
#Wilcoxon-Mann-Whitney秩和检验
wilcox_test(Prob~So,
data=UScrime,distribution="exact")
#范例三 k样本检验问题
library(mvtnorm)
library(survival)
library(TH.data)
library(MASS)
library(multcomp)
library(coin)
set.seed(1234)
#参考来自数据9999次的置换
#设定随机数种子可让结果重现
#结果表明各组间病人的响应值显著不同
oneway_test(response~trt,data=cholesterol,
distribution=approximate(B=9999))
#factor
f1 <- c("Type1","Type2","Type3","Type4")
f1 <- as.factor(f1)
class(f1) # 返回"factor"
f2 <- c(1,2,3,4)
f2 <- as.factor(f2)
class(f2) #返回 "factor"
#列联表中的独立性
library(coin)
library(vcd)
class(Arthritis$Improved) #返回 "ordered" "factor"
#为什么需要把变量Improved从一个有序因子变成一个分类因子
#因为,如果你用有序因子,coin将会生成一个线性与线性趋势检验,而不是卡方检验
#Approximative Linear-by-Linear Association Test
Ar<- transform(Arthritis,
Improved=as.factor(Improved))
class(Ar$Improved) #返回 "ordered" "factor" ,直接转换失败
Ar$Improved <- factor(Ar$Improved,ordered=FALSE) #必须加ordered=FALSE,才能转换成功
class(Ar$Improved) #返回 "factor" ,无序因子
set.seed(1234)
#Pearson 卡方检验
#Approximative Pearson Chi-Squared Test
chisq_test(Treatment~Improved,data=Ar,distribution=approximate(B=9999))
#p-value = 0.0006001,落入拒绝域,表明存在显著性差异
#数值变量间的独立性
#spearman_test()函数提供了两数值变量的独立性置换检验
#范例一
#state.x77原始是一个矩阵,在spearman_test中必须被转换成数据框
states <- as.data.frame(state.x77)
set.seed(1234)
spearman_test(Illiteracy~Murder,data=states,
distribution=approximate(B=9999))
#基于9999次重复的近似置换检验可知,独立性假设并不被满足
#两样本和k样本相关性检验
#两配对组的置换检验,可使用wilcoxsign_test()函数
#多于两组时,使用friedman_test函数
#范例一
library(coin)
library(MASS)
UScrime$U1 #14~24年龄段的失业率
UScrime$U2 #35~39年龄段的失业率
wilcoxsign_test(U1~U2,data=UScrime,distribution="exact")
#p-value = 1.421e-14,落入拒绝域,组别之间存在显著性差异
#lmPerm包
#专门用来做方差分析和回归分析的置换检验
#安装包
#下载lmPerm_1.1-2.tar.gz文件,下载路径:https://cran.r-project.org/src/contrib/Archive/lmPerm/
#安装Rtools 下载路径:https://cran.r-project.org/bin/windows/Rtools/Rtools35.exe
#install.packages(path,repos=NULL,type="source")
install.packages("D:\\Script\\R\\package\\lmPerm_1.1-2.tar.gz",repos=NULL,type="source")
#简单回归和多项式回归
#范例一
library(lmPerm)
set.seed(1234)
#lmp()即lm()的修改版,参数与lm()类似,仅额外多了一个perm参数:
#Exact:根据所有可能的排列组合生成精确检验
#Prob:从所有可能的排列中不断抽样,知道估计的标准差在估计的p值0.1以下,判停规则由可选的Ca参数控制
#若观测数大于10,perm="Exact",自动默认转为perm="Prob"
fit <- lmp(weight~height,data=women,perm="Prob")
summary(fit)
#增添的Iter栏中列出了要达到判停准则所需的迭代次数
#范例二
library(lmPerm)
set.seed(1234)
fit <- lmp(weight ~ height + I(height^2),data=women,perm="Prob")
summary(fit)
#多元回归
library(lmPerm)
set.seed(1234)
states <- as.data.frame(state.x77)
fit <- lmp(Murder~Population+Illiteracy+Income+Frost,
data=states,perm="Prob")
summary(fit)
#第8章中,正态理论中Population和Illiteracy均显著(P<0.05).
#而该置换检验中,Population不再显著
#当两种方法所得结果不一致时,你需要更加谨慎的审视数据,这很可能因为违反了正态性假设,或者存在离群点
#单因素方差分析
#范例一
#aovp()即avo()的修改版,参数与aov()类似,仅额外多了一个perm参数:
#Exact:根据所有可能的排列组合生成精确检验
#Prob:从所有可能的排列中不断抽样,知道估计的标准差在估计的p值0.1以下,判停规则由可选的Ca参数控制
#若观测数大于10,perm="Exact",自动默认转为perm="Prob"
library(lmPerm)
library(multcomp)
set.seed(1234)
fit <- aovp(response~trt,data=cholesterol,perm="Prob")
summary(fit)
anova(fit)
#Pr(Prob) < 2.2e-16 *** ,落入拒绝域,表明各疗法的效果不全相同
#范例二
library(lmPerm)
set.seed(1234)
fit <- aovp(weight~gesttime+dose,data=litter,perm="Prob")
anova(fit)
#依据p值,当控制gesttime,四种药物剂量对鼠患的体重影响不相同
#双因素方差分析
#当aovp应用到方差分析设计时,它默认使用唯一平方和方,也称为SAS三型平方和
#每种效应都会依据其他效应做响应调整。
#R中默认的参数化方差分析设计使用的是序贯平方和,SAS一型平方和
#对于平衡设计,两种方法结果相同
#对于不平衡设计,两种方法结果不相同,不平衡性越大,结果分歧越大
library(lmPerm)
set.seed(1234)
fit <- aovp(len~supp*dose,data=ToothGrowth,perm="Prob")
anova(fit)
#在0.05的显著性水平下,三种效应都不等于0;
#在0.01的水平下,只有主效应显著
#自助法
#所谓自助法,即从初始样本重复随机替换样本,生成一个或一系列待检验统计量的经验分布
#无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设
#自助法步骤
#1从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一直都不会被选中
#2计算并记录样本均值
#3重复1-2,1000次
#4将1000个样本均值从小到大排序
#5找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限定了95%的置信区间
#对单个统计量使用自助法
#范例一
#自助法的三个主要步骤
#步骤一:
#写一个能返回待研究统计量值的函数。
#如果只有单个统计量(如中位数),函数应该返回一个数值;
#如果有一列统计量(如一列回归系数),函数应该返回一个向量
#自定义获取R方的函数
rsq <- function(formula,data,indices){
d <- data[indices,] #必须申明,因为boot()要用其来选择样本
fit <- lm(formula,data=d)
return(summary(fit)$r.square)
}
#步骤二:
#为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理
library(boot)
set.seed(1234)
#boot(data=,statistic=,R=...)
#boot函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指标,
#这些指标被统计量函数用来选择样本
#data 向量,矩阵或者数据框
#statistic 生成k个统计量以供自举的函数,k=1时对单个统计量进行自助抽样
#R 自助抽样的次数
#boot函数对于rsq,相当调用rsq,1000次,rsq的形参data=mtcars,formula=mpg~wt+disp,indices=R=1000
results <- boot(data=mtcars,statistic=rsq,
R=1000,formula=mpg~wt+disp)
class(results) #返回"boot"
results$t0 #从原始数据得到的k个统计量的观测值,R方=0.7809306,等价于summary(fit)$r.square
results$t #一个R*K的矩阵,每行即k个统计量的自助重复值,1000个自助R方
#
print(results)
plot(results)
#步骤三:
#boot.ci(bootobject,conf=,type=)函数获取统计量的置信区间
#conf 预期的置信区间,默认conf=0.95
#type 获取置信区间的方法。可能值:
#norm,basic,stud,
#perc (分位数展示的是样本均值),
#bca (将根据偏差对区间做简单调整)
#all
boot.ci(results,type=c("perc","bca"))
#perc 对应 结果集Percentile数据列,代表用perc方法得到的95%置信区间的值
#bca 对应 结果集BCa数据列,代表用BCa方法得到的95%置信区间的值
#结论
#生成置信区间的不同方法将会导致获得不同的区间
#由于0都在置信区间以为,零假设H0,R平方值=0都被拒绝
#多个统计量的自助法
#范例一
#步骤一:
#自定义返回待研究统计量值的函数
bs <- function(formula,data,indices){
d <- data[indices,]
fit <- lm(formula,data=d)
return(coef(fit))
}#返回线性回归模型的回归系数
#步骤二
#通过boot达到对bs函数重复调用
library(boot)
set.seed(1234)
results <- boot(data=mtcars,statistic=bs,
R=1000,formula=mpg~wt+disp)
results$t0
#返回coef(fit)的实际值
#Intercept=34.96055
#wt=-3.350825
#disp=-0.01772474
#一个1000*3的矩阵,每行即3个统计量的自助重复值,1000个自助R方
results$t
#等价于results命令直接运行
#original列,等同于第一步自定义函数的真实统计量返回值
print(results)
plot(results,index=2) #绘制wt的直方图,散点图
#步骤三
#boot.ci提取置信区间
fit <- lm(mpg~wt+disp,data=mtcars)
class(coef(fit)) #返回numeric向量
coef(fit)[1] #返回截距回归系数 Intercept=34.96055
coef(fit)[2] #返回wt回归系数 wt=-3.350825
coef(fit)[3] #返回disp回归系数 disp=-0.01772474
summary(fit)$r.square #返回 r.square=0.7809306
boot.ci(results,type="bca",index=2) #利用bca方法,获取wt的95%的置信区间
boot.ci(results,type="bca",index=3) #利用bca方法,获取disp的95%的置信区间
#结论:
#一系列基于随机化和重抽样的计算机密集型方法,它们使你无需理论分布的知识便能够进行假设检验,
#获得置信区间。
#当数据来自于未知分布,或者存在严重的离群点,或者样本量过小,又或者没有参数方法可以回答你感兴趣的
#假设问题时,这些方法是非常实用的