读书笔记_第十二章

#重抽样与自助法

#置换检验
#又称为随机化检验，或重随机化检验
#假如你觉的假定数据成正态分布并不合适，或者担心离群点的影响，又或者感觉对于
#标准的参数方法来说数据集太小，那么置换检验便提供了一个非常不错的选择
#置换检验主要用于生成检验零假设的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%的置信区间

#结论：
#一系列基于随机化和重抽样的计算机密集型方法，它们使你无需理论分布的知识便能够进行假设检验，
#获得置信区间。
#当数据来自于未知分布，或者存在严重的离群点，或者样本量过小，又或者没有参数方法可以回答你感兴趣的
#假设问题时，这些方法是非常实用的