zoukankan      html  css  js  c++  java
  • 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)
    
    #置换检验和自助法不是万能的,烂数据无法转化为好数据。
    
  • 相关阅读:
    Linux命令应用大词典-第11章 Shell编程
    Kubernetes 学习12 kubernetes 存储卷
    linux dd命令
    Kubernetes 学习11 kubernetes ingress及ingress controller
    Kubernetes 学习10 Service资源
    Kubernetes 学习9 Pod控制器
    Kubernetes 学习8 Pod控制器
    Kubernetes 学习7 Pod控制器应用进阶2
    Kubernetes 学习6 Pod控制器应用进阶
    Kubernetes 学习5 kubernetes资源清单定义入门
  • 原文地址:https://www.cnblogs.com/jessepeng/p/10657603.html
Copyright © 2011-2022 走看看