R 语言实战-Part 3 笔记
R 语言实战(第二版)
part 3 中级方法
-------------第8章 回归------------------
#概念:用一个或多个自变量(预测变量)来预测因变量(响应变量)的方法
#最常用:OLS——普通最小二乘回归法,包括简单线性回归、多项式回归、多元线性回归
#过程:拟合OLS回归模型——>评价拟合优度——>假设检验——>选择模型
#OLS回归
#目标:减少因变量的真实值和预测值的差值来获得模型参数(截距和斜率),即使得残差平方和最小
#数据需满足:正态性、独立性、线性、同方差性
myfit <- lm(formula,data)
#表达式符号:~/+/:/*/^/./-1/I()/function
#拟合线性模型相关函数:
summary() #拟合模型详细结果
coefficients() #模型参数(截距和斜率)
confint() #模型参数置信区间(默认95%)
fitted() #模型预测值
residuals() #模型残差值
anova() #方差分析表
vcov() #模型参数的协方差矩阵
plot() #评价拟合模型诊断图
predict() #对新的数据集预测因变量
##1.简单线性回归
head(women)
fit <- lm(weight~height,data=women)
summary(fit) #理解结果
women$weight
fitted(fit)
residuals(fit) #真实值与预测值之差
plot(women$height,women$weight)
abline(fit)
#模型等式:weight=-87.52+3.45*height
#图形表明可用含一个弯曲的曲线来提高预测精度
##2.多项式回归
#注意不同于多元线性,它是n阶多项式,如x^2,x^3等,或者说它是多元线性回归的一个特例。
#多项式等式仍可认为是线性回归模型(加权和形式),非线性模型可用nls()函数拟合
#上例可通过添加一个二次项(即x^2)来提高回归模型精度
fit2 <- lm(weight~height+I(height^2),data = women)
summary(fit2)
plot(women$height,women$weight)
lines(women$height,fitted(fit2))
#p值、方差解释率都比原来模型要好。模型等式:weight=261.88-7.35*height+0.083*height^2
fit3 <- lm(weight~height+I(height^2)+I(height^3),data = women)
plot(women$height,women$weight)
lines(women$height,fitted(fit3))
#一般而言,n次多项式生成有n-1个弯曲的曲线,比三次更高的多项式几乎没有必要。
#快速绘制二元关系图
library(car)
scatterplot(weight~height,data=women,
spread=F, #删除残差正负均方根在平滑曲线的展开和对称信息
smoother.args=list(lty=2)) #设置拟合曲线为虚线
##3.多元线性回归
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
head(states)
##step 1:检查变量间相关性
cor(states)
car::scatterplotMatrix(states,spread=F,smoother.args=list(lty=2)) #因变量和自变量的散点图矩阵
##step 2:拟合多元线性回归模型
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
summary(fit)
#自变量多个时,回归系数含义为:一个自变量增加一个单位,其他自变量保持不变时,因变量将要增加的数量
##有交互项的多元线性回归
fit <- lm(mpg~hp+wt+hp:wt,data = mtcars) #马力与车重
summary(fit)
#两个自变量交互项显著,说明因变量与其中一个自变量的关系依赖于另一个自变量
#交互项展示: 看二者关系到底如何变化
library(effects)
plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline = T)
##4.回归诊断
#对模型参数推断的信心依赖于在多大程度上满足OLS模型统计假设
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
confint(fit)
#文盲率改变1%,谋杀率在95%的置信区间[2,38,5,90]内变化,Frost置信区间包含0,温度与谋杀率无关,但这些必须基于数据满足统计假设的前提。
##1)标准方法
#最常用:对lm函数返回对象plot得到的四幅图,这四幅图反应了OLS回归的统计假设:正态性、独立性、线性、同方差性
fit <- lm(weight~height,data=women)
par(mfrow=c(2,2))
plot(fit)
#正态Q-Q图:正态性,需落在45度角直线上
#残差图和拟合图residuals vs fitted:线性,需无关联
#位置尺度图scale-location graph: 同方差性,点应随机分布
#残差与杠杆图residuals vs leverage: 离群点、杠杆值、强影响点
##再看看二次拟合诊断图:
fit2 <- lm(weight~height+I(height^2),data=women)
par(mfrow=c(2,2))
plot(fit2)
#比之前理想一些
##剔除异常值(需谨慎)
newfit <- lm(weight~height+I(height^2),data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
par(mfrow=c(2,2))
plot(fit)
#剔除Nevada离群点,模型假设更好的满足
##2)改进方法
#car包titig提供了大量回归诊断函数
library(car)
qqPlot #分位数比较图
#正态性检验
qqPlot(fit,labels=row.names(states),
id.method="identify", #交互式绘图,用鼠标点,如实时查看异常值。按ESC,stop或右击图形退出该模式
simulate = T,main = "Q-Q plot")
#误差的独立性检验
durbinWatsonTest(fit) #若p不显著,则误差项之间独立
#自变量和因变量的线性检验
crPlots(fit) #成分残差图
#若图形不是线性,说明需要添加一些曲线信息,如多项式或log转换等
#同方差性检验
ncvTest(fit) #若p不显著,说明方差不变
spreadLevelPlot(fit) #若方差稳定,将会一些水平趋势线
##线性模型的综合检验
library(gvlma)
summary(gvlma(fit))
##3)异常观测值
###离群点(预测不佳的观测点)
outlierTest(fit) #根据p值一次只能返回最大的那个离群点,只能删了该值才能继续检测
###高杠杆值点(许多异常预测变量的组合)
#帽子统计量来判断
###强影响点(对模型参数估计比例失衡的点)
#通过Cook距离(D统计量)或变量添加图来检测
###整合离群点、杠杆值和强影响点到一幅图中
influencePlot(fit)
##5.改进模型
#删除观测点(谨慎):重新拟合
#变量变换:不符合正态性、线性或同方差性假设时。PS:若变换了变量,解释必须基于变换后的变量,而非初始变量
#增删变量:处理多重共线性时
#尝试其他方法:异常值——稳健回归模型,非正态——非参数回归模型,非线性——非线性回归模型,不同方差——时间序列模型或多层次回归模型
##6.选择最佳的模型
#回归模型的选择总会涉及到:预测精度和模型简洁度的调和问题
##1)两个模型比较
#anova方法:需要嵌套模型
fit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
fit2 <- lm(Murder~Population,data = states) #fit1和fit2是嵌套模型(一个的项包含另一个)
anova(fit1,fit2) #不显著,说明可以将后两项不加入线性模型中
#AIC方法:不需嵌套模型
AIC(fit1,fit2) #直接根据AIC值判断(越小越好)
##2)变量选择
#逐步回归法:
#模型一次添加(向前)或删除(向后)一个变量,直至判停。通常结合两者。
#可能会找到一个好模型,但不能保证是最佳模型,因为没有评价所有可能模型。
library(MASS)
stepAIC(fit,direction = "backward") #变量从多到少,直到找到最小的AIC
#全子集回归:
#所有可能的模型都会被检验。R平方含义是自变量解释因变量的程度,调整R平方考虑了模型参数数目。
#一般要优于逐步回归,但预测变量太多时,速度很慢。
#方法1
library(leaps)
leaps <- regsubsets(Murder~Population+Illiteracy+Income+Frost,data = states,nbest = 4)
plot(leaps,scale = "adjr2") #纵轴调整R平方
#方法2
library(car)
subsets(leaps,statistic = "cp") #交互绘图,鼠标点图例
abline(1,1,lty=2,col="red") #模型越好,离这条直线越近
##6.交叉验证
#评价回归方程的泛化能力
#概念:挑出一定比例数据作为训练样本,其他样本作保留样本,先在训练样本上获得回归方程,再在保留样本上预测。
#k重交叉验证:样本被分为k个子集,轮流将k-1个子集组合作为训练集,另一个子集作为保留集,这样得到的k个预测方程,记录k个保留集的预测结果,再求均值。
#bootstrap::crossval函数
shrinkage <- function(fit,k=10){
library(bootstrap)
theta.fit <- function(x,y){lsfit(x,y)}
theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}
x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]
results <- crossval(x,y,theta.fit,theta.predict,ngroup=k)
r2 <- cor(y,fit$fitted.values)^2
r2cv <- cor(y,results$cv.fit)^2
cat("original R square = ",r2,"\n")
cat(k,"fold cross validated R square =",r2cv,"\n")
cat("change =",r2-r2cv,"\n")
}
shrinkage(fit1) #每次都会有少许不同
#比较验证前R平方和交叉验证后R平方,R平方减少得越少,预测越精确
shrinkage(fit2)
#fit2模型的R平方减少的更少,可见fit2模型更好
##变量的相对重要性
#大部分情况下,自变量(预测变量)之间都有一定相关性,因此需要对自变量的重要性进行排序,即比较标准化后的回归系数
zstates <- as.data.frame(scale(states)) #lm函数需要数据框格式
zfit <- lm(Murder~Population+Illiteracy+Income+Frost,data = zstates)
coef(zfit) #可以看出,Illiteracy最大
##相对权重方法(relwights函数代码见书)
relweights(fit)
##总结:回归分析是很多方法的总称,包括拟合模型、假设检验、修正模型、再拟合等过程
---------------------第9章 方差分析---------------------------------------------
##1)术语
#anova
#单因素方差分析:1个类别型变量,分为组间(因素)和组内(水平),通过F检验评价,均衡设计(观测数相同)和非均衡设计
#双因素方差分析:两个因子(如时间和疗法),主效应(两个因子的影响),交互效应(交互影响),双因素混合模型方差分析(包括组间和组内设计)
#多因素方差分析:多个因子
#协方差分析:ancova,即存在混淆因素,干扰变数,如抑郁症干扰焦虑症
#多元方差分析:manova,2个及其以上个因变量
#多元协方差分析:mancova,即存在协变量的多元方差分析
##2)anova模型拟合
#方差分析是广义线性模型的特例,lm函数也能分析,aov函数与之等同
#表达式:小写字母定量,大写字母组别
y~A #单因素anova
y~x+A #单个协变量的单因素ancova
y~A*B #双因素anova
y~x1+x2+A*B #含两个协变量的 双因素ancova
y~B+A #随机化区组(B是区组因子)
y~A+Error(Subject/A) #单因素组内anova
y~B*W+Error(Subject/W) #含单个组内因子(W)和单个组间因子(B)的重复测量anova
#R默认序贯型方法计算anova效应:y~A+B+A:B(A对y的影响;控制A时,B对y的影响;控制A和B主效应时,A和B的交互效应)
#样本大小越不平衡,效应项的顺序对结果影响越大
#基本准则:若研究设计不是正交的(因子和协变量相关),一定要谨慎设置效应的顺序
##3)单因素方差分析
library(multcomp)
attach(cholesterol)
cholesterol #五种疗法
table(trt)
aggregate(response,by=list(trt),mean)
aggregate(response,by=list(trt),sd)
fit <- aov(response~trt)
summary(fit) #p<0.01,F检验显著,说明五种疗法效果不同,但不知哪种显著差异
library(gplots)
plotmeans(response~trt) #各组均值及其置信区间
detach(cholesterol)
##多重比较
###方法1:
TukeyHSD(fit) #各组比较结果
par(las=2) #旋转轴标签
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit)) #置信区间包含0的说明不显著
###方法2:
tuk <- multcomp::glht(fit,linfct=mcp(trt="Tukey"))
plot(cld(tuk,level = 0.05), #0.05显著水平(95%置信区间)
col="lightgrey") #有相同字母的均值差异不显著
##方差分析假设检验条件
###正态性:qq图
qqPlot(lm(response~trt,data = cholesterol),
simulate=T) #数据落在95%置信区间范围,满足正态性(只能大致判断)
###方差齐性:
bartlett.test(response~trt,data = cholesterol) #各组方差无显著性差异,满足条件
#方差齐性对离群点很敏感,检测离群点:
car::outlierTest(fit) #p>1时产生NA
##3)单因素协方差分析
data(litter,package = "multcomp")#4组怀孕小鼠,用4种不同剂量的药物对幼崽出生体重的影响,怀孕时间为协变量
help(litter,package = "multcomp")
attach(litter)
table(dose)
aggregate(weight,by=list(dose),mean)
fit <- aov(weight~gesttime+dose)
summary(fit)
#F检验怀孕时间和药物剂量都显著
#去除协变量效应后的组均值
library(effects)
effect("dose",fit)
#多重比较
multcomp::glht()
#假设检验条件
#方差齐性和正态性检验与方差分析类似,另外还要检验回归斜率的同质性
fit2 <- aov(weight~gesttime*dose,data = litter)
summary(fit2)
#交互效应不显著,说明斜率相等
#可视化
library(HH)
ancova(weight~gesttime+dose,data=litter)
#回归线(斜率)相互平行,截距不同
##4)双因素方差分析
help("ToothGrowth")
head(ToothGrowth)#两种喂食方法,分别3种抗坏血酸含量水平,对牙齿长度影响
attach(ToothGrowth)
table(supp,dose)
aggregate(len,by=list(supp,dose),mean)
aggregate(len,by=list(supp,dose),sd)
fit <- aov(len~supp*as.factor(dose)) #剂量要转换为因子变量,这样才是分组变量,而非数值型协变量
summary(fit) #主效应和交互效应都显著
#交互效应可视化
interaction.plot(dose,supp,len,type = "b",col = c("red","blue"),pch = c(16,18))
#或
gplots::plotmeans(len~interaction(supp,dose,sep=" "),
connect = list(c(1,3,5),c(2,4,6)),
col = c("red","darkgreen")) #展示了均值、误差棒(95%置信区间)和样本大小
#所有结果可视化(包括主效应和交互效应)
HH::interaction2wt(len~supp*dose) #同样可用于展示多因素方差分析
detach(ToothGrowth)
##5)重复测量方差分析、多元方差分析 略
##6)用回归来做anova
#anova和回归都是广义线性模型的特例
library(multcomp)
head(cholesterol)
levels(cholesterol$trt)
#anova拟合
summary(aov(response~trt,data = cholesterol))
#lm拟合
summary(lm(response~trt,data = cholesterol))
##线性模型要求自变量为数值型,所以当lm碰到因子时,会用与因子水平对应的数值型来代替因子,如trt有5个因子水平,它会创建4个对照变量
##这里对照没搞懂?
-------------------------第10章 功效分析-----------------------------
#即在一定置信度下,需要多少样本量才能检测到效应值?
#功效分析针对的是假设检验
##1)假设检验基本概念
#零假设H0,备择假设H1,显著性水平。样本统计量来推断总体参数
#I型错误(拒绝了真实的零假设,假阳性错误),II型错误(没拒绝错误的零假设,假阴性错误)
#关注4个量:给定任意3个,可推算第4个
#样本大小——观测数目
#显著性水平——假阳性概率
#功效——1-假阴性概率
#效应值——备择假设下效应的量
##2)pwr包做功效分析
library(pwr)
#一般先根据各种检验的各自方法计算出效应值,再计算所需样本量
#①t检验
pwr.t.test(n = ,d = ,sig.level = ,power = ,type = ,alternative = )
##n样本大小,d为功效值(标准化的均值之差),sig.level显著性水平,power功效水平,type检验类型(单样本/双样本),alternative双侧/单侧检验(less/greater)
pwr.t.test(d=.8,sig.level = .05,power = .9,type = "two.sample",alternative = "two.sided")
#计算n样本(每组需要的样本量)
pwr.t.test(n = 20,d=.5,sig.level = .01,type = "two.sample",alternative = "two.sided")
#计算功效水平power
#两样本大小不等的t检验
pwr.t2n.test(n1=,n2=,d=,sig.level = ,power = ,alternative = )
#②方差分析
pwr.anova.test(k=,n=,f=,sig.level = ,power = )
#k是组的个数,n是各组样本大小,f是效应值(公式略),sig.level显著水平,power功效水平
pwr.anova.test(k=5,f=0.25,sig.level = 0.05,power = 0.8)
#计算各组需要的样本量n
#③相关性
pwr.r.test(n=,r=,sig.level = ,power = ,alternative = )
#r为效应值(线性相关系数)
pwr.r.test(r=0.25,sig.level = 0.05,power = 0.9,alternative = "greater")
#计算样本量n
#④线性模型
pwr.f2.test(u=,v=,f2=,sig.level = ,power = )
#u分子自由度,v分母自由度,f2是效应值(公式有两种,略)
pwr.f2.test(u=3,f2=0.0769,sig.level = 0.05,power = 0.9)
#计算样本量
#⑤比例检验
pwr.2p.test(h=,n=,sig.level = ,power = ) #h为效应值,n为各组相同的样本量
#⑥卡方检验:一般用来评价两个类别型变量的关系(是否独立)
pwr.chisq.test(w=,N=,df=,sig.level = ,power = )
#w效应值,N总样本量,df自由度
#功效分析中,预期效应值是最难决定的参数,需要对主题有一定的经验和了解。
#Cohen效应值基准:可划分为小、中、大三种效应值
#⑦绘制功效分析图形
#检验各种效应值下相关性所需的样本量曲线
library(pwr)
#生成一系列相关系数和功效值
r <- seq(0.1,0.5,0.01) #相关系数
nr <- length(r)
p <- seq(0.4,0.9,0.1) #功效值
np <- length(p)
#获取样本大小
samsize <- array(numeric(nr*np),dim = c(nr,np)) #先建立空的样本矩阵,再用循环填充
for(i in 1:np){
for (j in 1:nr) {
result <- pwr.r.test(n=NULL,r=r[j],sig.level = 0.05,power = p[i],alternative = "two.sided")
samsize[j,i] <- ceiling(result$n) #向上取整
}
}
#创建图形
xrange <- range(r)
yrange <- round(range(samsize))
colors <- rainbow(length(p))
plot(xrange,yrange,type = "n")
#添加功效曲线
for (i in 1:np) {
lines(r,samsize[,i],type = "l",lwd=2,col=colors[i])
}
#添加网格线
abline(v=0,h=seq(0,yrange[2],50),lty=2,col="grey89")
abline(h=0,v=seq(xrange[1],xrange[2],0.02),lty=2,col="gray89")
#添加注释
title("test\nSig=0.05 (Two-tailed)")
legend("topright",title = "power",as.character(p),fill = colors)
#基因研究中功效分析的R包:
#powerGWASinteraction pedantics gap ssize.fdr
----------------------第11章 中级绘图--------------------------------
##二元变量和多元变量关系的可视化
#1.散点图
##基本散点图
attach(mtcars)
plot(wt,mpg,
main = "scater plot of mpg vs weight",
xlab = "car weight",
ylab = "miles per gallon",
pch=19
)
abline(lm(mpg~wt),col="red",lwd=2,lty=1) #添加线性拟合直线
lines(lowess(wt,mpg),col="blue",lwd=2,lty=2) #添加lowess拟合曲线(R语言中还有一个loess平滑曲线拟合函数)
library(car)
scatterplot(mpg~wt | cyl,data=mtcars, #按条件绘图(即按cyl的水平分别绘制前两者关系图)
lwd=2,span=0.75, #span控制loess曲线平滑量
legend.plot=T,
id.methods="identify", #鼠标单击交互
#labels=row.names(mtcars),
boxplots="xy") #边界箱线图
##散点图矩阵
pairs(~mpg+disp+drat+wt,data=mtcars)
pairs(~mpg+disp+drat+wt,data=mtcars,upper.panel=NULL) #只展示一边(另一半完全相同)
pairs(~mpg+disp+drat+wt,data=mtcars,lower.panel=NULL)
library(car)
scatterplotMatrix(~mpg+disp+drat+wt,data=mtcars,
spread=F, #不添加展示分散度和对称信息的直线
smoother.args=list(lty=2)) #平滑(loess)拟合曲线为虚线
#为何每个图有三条平滑曲线?
##高密度散点图
#当数据点重叠很严重时用普通散点图很难识别哪里数据点密度大
set.seed(123)
mydata <- as.data.frame(rbind(matrix(rnorm(10000,0,0.5),ncol = 2),
matrix(rnorm(10000,3,2),ncol = 2)))
names(mydata) <- c("x","y")
head(mydata)
with(mydata,plot(x,y,pch=19))
#smoothScatter用核密度估计生成颜色密度
with(mydata,smoothScatter(x,y))
library(hexbin)
with(mydata,plot(hexbin(x,y,xbins=50)))
##三维散点图
#三个定量变量的交互关系
attach(mtcars)
library(scatterplot3d)
s3d <- scatterplot3d(wt,disp,mpg,
pch=16,
highlight.3d = T,
type = "h")
s3d$plane3d(lm(mpg~wt+disp)) #添加一个回归面(代表预测值)
##旋转三维散点图: 可交互
library(rgl)
plot3d(wt,disp,mpg,col="red",size=5)
library(car)
with(mtcars,scatter3d(wt,disp,mpg))
##气泡图
#创建二维散点图,再用点大小表示第三个变量
symbols(wt,mpg,circle=sqrt(disp/pi),
inches = 0.3, #比例因子
fg="white",bg="lightblue")
text(wt,mpg,rownames(mtcars),cex = 0.6)
#2.折线图
##两个函数plot(x,y,type=)或lines(x,y,type=)
help(Orange)
opar <- par(no.readonly = T)
par(mfrow=c(1,2))
t1 <- subset(Orange,Tree==1)
plot(t1$age,t1$circumference)
plot(t1$age,t1$circumference,type="b") #类型有p/l/o/b/n/c/s/h等
par(opar)
#可用plot函数先创建空图形(轴标签、轴范围等),再用lines函数添加数据点
#例子:展示五种橘树随时间推移的生长状况的折线图 <学习编程思维>
Orange$Tree <- as.numeric(Orange$Tree)
ntrees <- max(Orange$Tree)
#创建图形
xrange <- range(Orange$age)
yrange <- range(Orange$circumference)
plot(xrange,yrange,type="n",
xlab="Age (days)",
ylab="Circumference (mm)")
colors <- rainbow(ntrees)
linetype <- c(1:ntrees)
plotchar <- seq(18,18+ntrees,1) #点类型
#添加线条
for (i in 1:ntrees) {
tree <- subset(Orange,Tree==i)
lines(tree$age,tree$circumference,
type = "b",lwd=2,
lty=linetype[i],col=colors[i],pch=plotchar[i])
}
title("Tree growth","example of line plot")
#添加图例
legend(xrange[1],yrange[2], #位置(但这里好像不对?)
1:ntrees,
cex = 0.8,title = "Tree",
col = colors,lty = linetype,pch = plotchar)
#3.相关图
options(digits = 2)
cor(mtcars)
library(corrgram)
corrgram(mtcars,order = T, #行列重排序,将相似相关的变量聚在一起(使用主成分法)
lower.panel = panel.shade, #阴影深度表相关性大小
upper.panel = panel.pie, #饼图比例表相关性大小
text.panel = panel.txt, #变量名
main="corrgram of mtcars inttercorrelations"
)
#蓝色(斜向上)正相关,红色(斜向下)负相关。还可对主对角线进行设置
corrgram(mtcars,order = T,lower.panel = panel.ellipse, #置信椭圆和平滑曲线
upper.panel = panel.pts, #散点图
diag.panel = panel.minmax, #输出最大最小值
text.panel = panel.txt)
corrgram(mtcars,lower.panel = panel.shade,
upper.panel = NULL, #只展示半边
text.panel = panel.txt)
corrgram(mtcars,lower.panel = panel.shade,order = T,
upper.panel = NULL, text.panel = panel.txt,
col.regions = colorRampPalette(c("red","yellow","blue","green"))) #指定颜色
#4.马赛克图
#展示多个类别型变量的关系
ftable(Titanic) #紧凑列联表
library(vcd)
mosaic(ftable(Titanic))
mosaic(Titanic,shade = T,legend=T)#泰坦尼克号船舱等级、性别、年龄和幸存者关系,信息量很大
#等同于
mosaic(~Class+Sex+Age+Survived,data=Titanic,shade = T,legend=T) #更可控个性化
mosaic(~Class+Sex+Survived,data=Titanic,shade = T,legend=T)
------------------------第12章 重抽样与自助法----------------------------------------
#当不满足统计假设(未知/混合分布、样本过小、离群点等)时,基于随机化和重抽样统计方法适用
#1.置换检验
#也叫随机化检验,将统计量与置换观测后的经验分布(而非理论分布)进行比较,根据统计量值的极端性判断是否拒绝零假设
#经验分布依据所有可能的排列组合,称为精确检验(只适合样本)
#从所有可能的排列中进行抽样获得近似的检验,称为蒙特卡洛模拟
#置换检验需要设定随机数种子,便于结果重现
#两个用于置换检验的R包:
#coin包对于独立性问题提供全面的置换检验的框架
#lmPerm包专门用来做方差分析和回归分析的置换检验
#2)coin包
#可解决独立性的问题:因变量和组别分配独立吗?两个数值变量独立吗?两个类别变量独立吗?
#coin包提供的函数与传统检验类似
##①两样本检验
library(coin)
score <- c(40,57,45,55,58,57,64,55,62,65)
treatment <- factor(c(rep("A",5),rep("B",5))) #coin包中类别和序数变量都要转化为因子
mydata <- data.frame(treatment,score)
#传统t检验
t.test(score~treatment,data = mydata,var.equal=T)
#置换检验
oneway_test(score~treatment,data = mydata,distribution="exact") #精确检验
#尽管t检验p值显著,置换检验不显著,但由于样本量小,反而更相信置换检验的结果
##②K样本置换检验
library(multcomp)
set.seed(1234)
oneway_test(response~trt,data = cholesterol,
distibution=approximate(B=9999)) #置换9999次,p显著
##③列联表的独立性
library(vcd)
Arthritis
Arthritis <- transform(Arthritis,Improved=as.factor(as.numeric(Improved)))
#将Improved从一个有序因子变成一个分类因子,因为是做卡方检验,而非趋势检验
set.seed(1234)
chisq_test(Treatment~Improved,data = Arthritis,
distribution=approximate(B=9999))
##④数值变量间的独立性
set.seed(123)
spearman_test(Illiteracy~Murder,data=as.data.frame(state.x77),
distribution=approximate(B=9999))
##⑤两样本和K样本相关性检验
#两配对样本置换检验:
library(MASS)
wilcoxsign_test(U1~U2,data=UScrime,distribution="exact")
#多于两组时:
friedman_test()
#3)lmPerm包
#用于做线性模型的置换检验
#lmp和aovp函数对应正态理论检验的lm和aov函数,多一个perm=参数,可选Exact(精确检验,样本量小于10)/Prob/SPR(贯序概率比)
##①简单回归和多项式回归
library(lmPerm)
set.seed(123)
fit <- lmp(weight~height,data = women,perm = "Prob") #简单回归
summary(fit)
set.seed(1234)
fit <- lmp(weight~height+I(height^2),data = women,perm = "Prob") #多项式回归
summary(fit)
#结果和lm函数类似,多一栏Iter列出要达到判停所需的迭代次数
##②多元回归
set.seed(12345)
fit <- lmp(Murder~Population+Illiteracy+Income+Frost,
data=as.data.frame(state.x77),perm = "Prob")
summary(fit)
#与第8章lm函数结果不同
##③单因素方差分析和协方差分析
library(multcomp)
set.seed(1234)
fit <- aovp(response~trt,data = cholesterol,perm = "Prob")
anova(fit)
fit <- aovp(weight~gesttime+dose,data = litter,perm = "Prob")
anova(fit)
##④双因素方差分析
fit <- aovp(len~supp*dose,data=ToothGrowth,perm = "Prob")
anova(fit)
#置换检验适合处理非正态数据、存在离群点、样本很小或者无法做参数检验等情况。有助于"效应是否存在"这样的问题。
#置换方法对于获取置信区间和估计测量精度比较困难,自助法派上用场。
#2.自助法
#初始样本重复随机抽样,生成一个或一系列待检验统计量的经验分布,无需假设理论分布,便可生成统计量的置信区间,并能检验统计假设
#boot包中的boot()函数和boot.ci()函数
#疑问:统计量有哪些?R平方,回归系数,summary函数结果等
##1)单个统计量的自助法
#获取R平方函数
rsq <- function(formula,data,indices){
fit <- lm(formula,data=data[indices,])
return(summary(fit)$r.square)
}
#做大量自助抽样
library(boot)
set.seed(123)
results <- boot(data=mtcars,statistic = rsq,
R=1000,formula=mpg~wt+disp) #自助1000次
#输出对象
print(results)
#绘制结果
plot(results)
#获取95%置信区间
boot.ci(results,type = c("perc","bca")) #perc和bca两种方法得到的结果略有不同
#2)多个统计量的自助法
bs <- function(formula,data,indices){
fit <- lm(formula,data=data[indices,])
return(coef(fit))
}
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)
boot.ci(results,type = "bca",index = 3)
#置换检验和自助法不是万能的,烂数据无法转化为好数据。
本文来自博客园,作者:生物信息与育种,转载请注明原文链接:https://www.cnblogs.com/miyuanbiotech/p/10657603.html。若要及时了解动态信息,请关注同名微信公众号:生物信息与育种。