R语言实战-第十二章重抽样与自助法

#第12章 重抽样与自助法
library(coin)
install.packages(file.choose(),repos = NULL,type = "source")
library(lmPerm)
#为什么选择置换检验 数据呈正态分布不太合适 或者担心离群点 
#又或者对标准的参数方法来说数据集太小
#当样本量较大时 可以使用蒙特卡洛模拟 进行抽样 获得一个近似的检验
#coin包对于独立性问题提供了一个很全面的置换检验框架
#lmPerm包专门用来做方差分析和回归分析的置换检验

#用coin包做置换检验 类别型变量和序数变量必须分别转化成因子和有序因子；数据以数据框形式储存

#两样本数相同 n1=n2=5
#独立两样本 t检验与单因素置换检验

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")
#t.test与oneway_test结果p值并不一致 由于数据量较小作者更相信置换检验结果 

#两样本数不同 n1=31 n2=16
#Wilcoxon-Mann-Whiney U检验 检验美国南部监禁概率与非南部间的差异

library(MASS)
with(UScrime,table(So))
UScrime <- transform(UScrime,So=factor(So))#因子转化
wilcox_test(Prob~So,data = UScrime,distribution="exact")

# 5组样本 样本数均为10
#K样本检验 评价五种药物疗法对降低胆固醇的效果

library(multcomp)
set.seed(1234)
with(cholesterol,table(trt))
oneway_test(response~trt,data = cholesterol,
            distribution=approximate(nresample=9999))

#列联表中的独立性

#判断两类别型变量的独立性 
#通过chisq_test()或cmh_test()函数
#若变量都是有序型 可以使用lbl_test()函数检验是否存在线性趋势
#检验治疗(2个水平)与效果(3个水平)的关系
library(coin)
library(vcd)
head(Arthritis)
Arthritis <- transform(Arthritis,Improved=as.factor(as.numeric(Improved)))#将效果转化成分类因子
str(Arthritis)
set.seed(1234)
chisq_test(Treatment~Improved,data = Arthritis,
           distribution=approximate(nresample=9999))
#若将Improved转化成有序因子而不是分类因子 coin()将会生成一个线性与线性趋势检验 而非卡方检验

#数值变量间的独立性
#两数值变量的独立性置换检验 spearman_test()函数
#检验美国文盲率与谋杀率之间的相关性
states <- as.data.frame(state.x77)#将矩阵转化成数据框
set.seed(1234)
spearman_test(Illiteracy~Murder,data = states,
              distribution=approximate(nresample=9999))

#两样本和K样本相关性检验
#对于两配对组的置换检验 使用wilcoxsign_test()函数
#多于两组时 使用friedman_test()函数
#两个年龄段间失业率差异
library(coin)
library(MASS)
library(Hmisc)
uuu <- c("U1","U2")
Hmisc::describe(UScrime[uuu])
wilcoxsign_test(U1~U2,data=UScrime,
                distribution="exact")
#结果表明两者失业率是不同的

#lmPerm包的置换检验 可以做线性模型的置换检验
#比如lmp()就是lm()的修改 aovp()就是aov()的修改版 能进行置换检验而非正态理论检验
#当观测数>10时 perm="exact"会自动默认转成perm="prob"
#简单回归和多项式回归

#简单线性回归的置换检验
#women数据集 分析体重与升高的相关性

library(lmPerm)
summary(women)
str(women)
set.seed(1234)
fit <- lmp(weight~height,data = women,perm = "Prob")
summary(fit)

#多项式回归的置换检验
#women数据集 分析体重与身高的相关性
set.seed(1234)
fit <- lmp(weight~height+I(height^2),data = women,perm = "Prob")
summary(fit)

#多元回归
#美国50个州的人口数 文盲率 收入水平 结霜天数预测犯罪率
states <- as.data.frame(state.x77)
set.seed(1234)
fit <- lmp(Murder~Population+Income+Illiteracy+Frost,
           data = states,perm = "Prob")
summary(fit)
#与正态理论中的结果不一样 这是因为存在离群点Nevada

#单因素方差分析和协方差分析
#单因素方差分析的置换检验
#疗法对降低胆固醇的影响
library(multcomp)
set.seed(1234)
fit <- aovp(response~trt,
           data = cholesterol,perm = "Prob")
summary(fit)
anova(fit)

#单因素协方差分析的置换检验
#当控制妊娠时间相同时 观测四种药物剂量对老鼠幼崽体重的影响
set.seed(1234)
fit <- aovp(weight~gesttime+dose,
            data = litter,perm = "Prob")
anova(fit)

#双因素方差分析
#喂食剂量和方式对豚鼠牙齿生长的影响
set.seed(1234)
fit <- aovp(len~supp*dose,
            data = ToothGrowth,perm = "Prob")
anova(fit)
#aovp()默认使用位移平方和法 即SAS类型III平方和

#置换检验点评
#置换检验为不理会正态分布 t分布 F分布或卡方分布提供了可能 
#可以处理非正态数据(比如分布偏倚很大) 存在离群点 样本量很小或无法做参数检验等情况
#置换方法无法获取置信区间和估计测量精度(可以使用自助法实现)
sqrt(10)
library(boot)
#对单个统计量使用自助法
#####预测变量对响应变量可解释的方差比
#写一个获取R平方的值的函数
rsq <- function(formula,data,indices){
  d <- data[indices,]
  fit <- lm(formula,data = d)
  return(summary(fit)$r.square)
}
#做大量的自助抽样 如1000次
library(boot)
set.seed(1234)
results <- boot(data=mtcars,statistic=rsq,
                R=1000,formula=mpg~wt+disp)
#boot的对象可以输出
print(results)
plot(results)

#获得95%的置信区间
boot.ci(results,type = c("perc","bca"))

#多个统计量的自助法
#####三个回归系数 截距 车重 发动机排量 95%的置信区间
#创建一个返回回归系数向量的函数
bs <- function(formula,data,indices){
  d <- data[indices,]
  fit <- lm(formula,data = d)
  return(coef(fit))
}

#使用该函数自助抽样1000次
library(boot)
set.seed(1234)
results <- boot(data=mtcars,statistic=bs,
                R=1000,formula=mpg~wt+disp)

print(results)
#index123分别是截距项 车重 发动机排量

#获得车重和发动机排量95%的置信区间
boot.ci(results,type = "bca",index = 2)#车重置信区间
boot.ci(results,type = "bca",index = 3)#发动机排量置信区间