zoukankan      html  css  js  c++  java
  • R语言--史上最全方差分析

    单因素

    rm(list = ls())
    install.packages('multcomp')
    library(multcomp)
    data('cholesterol')
    head(cholesterol)
    data <- cholesterol
    #正态性分析####
    shapiro.test(cholesterol$response) #data#观测值
    #不满足正态性
    #1.当偏态分布时,且偏态分布不严重时,仍可以使用方差分析,应当报告偏度、峰度
    fBasics::basicStats(x)
    #2.使用非参数检验Kruskal-Waliis检验
    kruskal.test()
    
    #方差齐性检验####
    bartlett.test(response~trt,data=cholesterol) #观测值~分组
    
    library(car)
    leveneTest(response ~ trt,data=cholesterol, center=mean)#观测值~分组
                                                #center=mean为原始levene,median更稳健
    
    #不满足方差齐性####
    #1.应先处理异常值
    #2.使用Welch's anova方法,两两比较使用Game-Howell方法(尤其各组例数不同)
    oneway.test(x~g,data=,var.equal = F)
    #Game-Howell方法:http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R
    
    #单因素方差分析
    
    anova <- with(cholesterol,aov(response ~ trt))
    summary(anova)
    
    OR
    
    anova <- oneway.test(response ~ trt , cholesterol, var.equal = T)
    anova
    ####描述统计
    by(cholesterol$response,cholesterol[,"trt"],summary) 
    aggregate(cholesterol$response,by=list(cholesterol$trt),FUN=sd)
    fBasics::basicStats(x)
    
    #自定义函数
    mystats <- function(x,na.omit=FALSE){   #一个函数返回x的多个描述性统计量
      if (na.omit)
        x <- x[!is.na(x)]
      m <- mean(x)
      n <- length(x)
      s <- sd(x)
      return(c(n=n,mean=m,stdev=s))
    }
    dstats <- function(x)sapply(x, mystats)
    
    options(digits = 3)
    by(data[2],data$trt,dstats)  
    #data$response 和 data[,2]数据均为横向显示,而data[2]则为纵向
    by(mtcars[1],mtcars$am,dstats)
    
    #画图,建议用graphpad prism画
    install.packages('gplots')
    library(gplots)
    attach(cholesterol)
    plotmeans(response ~ trt)
    detach(cholesterol)
    
    ##单因素方差分析与单因素线性回归的关系
    anolm <-lm(response ~trt ,data = cholesterol)  #trt为有4个哑变量,把1time当做reference
    summary(anolm) #F值相等

    两因素

    rm(list = ls())
    data("ToothGrowth")
    head(ToothGrowth)
    
    table(ToothGrowth$supp,ToothGrowth$dose)  #查看各组例数
    
    #各组均数 标准差
    aggregate(ToothGrowth$len,by=list(ToothGrowth$supp,ToothGrowth$dose),mean)
    aggregate(ToothGrowth$len,by=list(ToothGrowth$supp,ToothGrowth$dose),sd)
    
    #双因素方差分析
    twoanova <- aov(len ~ supp*dose,data = ToothGrowth)#+只有主效应| :只有交互效应|*都有
    summary(twoanova)
    
    #画图
    interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,
                     type='o',
                     col=c('red','blue'))

    两两比较

    rm(list = ls())
    
    data("airquality")
    head(airquality)
    
    #两两比较
    #两两比较方法很多,具体的选择看个人
    #最后的个人小结(仅代表个人观点): 如果各组例数相等,建议首选Tukey法;
    #如果例数不等,建议首选Scheffe法(如果比较组数不多,如3组,Bonferroni法也可以作为首选);
    #如果要分别比较每个试验组与对照组,建议采用Dunnett法;
    
    model <- aov(Temp~as.factor(Month),data=airquality)
    summary(model)
    ##Tukey
    
    TukeyHSD(model)
    plot(TukeyHSD(model))#除了横跨虚线的均有统计学意义
    
    
    ##LSD医学常用
    install.packages("agricolae")
    library(agricolae)
    
    out <- LSD.test(y = model, trt = "as.factor(Month)", p.adj = "none" ) 
                                                #多重比较,不矫正P值,==t检验
                                                #p.adj="bonferroni"
    out$groups # 查看每个组的label,group不一样即有差别,
    
    plot(out) # 可视化展示
    
    
    ##如果例数不等,建议首选Scheffe法
    library(agricolae)
    
    out <-scheffe.test (y = model,trt = "as.factor(Month)")
    
    out$group 
    
    plot(out) 
    
    
    ##如果要分别比较每个试验组与对照组,建议采用Dunnett法
    library(multcomp)
    
    out <- glht(model1, 
                linfct = mcp('as.factor(Month)' = 'Dunnett'), #mcp可以是不同的方法
                alternative = 'two.side')
    
    summary(out)

    协方差

    rm(list = ls())
    library(multcomp)
    data(litter)
    head(litter)
    ###协方差分析指的是某一变量A影响了自变量对响应变量的效应,变量A被称作协变量
    #前提要求:协变量为连续性变量 | 响应变量正态且方差齐|协变量对响应变量的影响呈线性关系
    
    
    #正态和方差齐
    table(litter$dose)
    aggregate(litter$weight,by=list(litter$dose),FUN=mean)
    aggregate(litter$weight,by=list(litter$dose),FUN=sd)
    shapiro.test(litter$weight)
    bartlett.test(litter$weight~litter$dose)
    library(car)
    leveneTest(weight ~ dose, data=litter)
    qqPlot(lm(weight ~ gesttime + dose,data=litter),simulate = T,labels=F)
    #协变量对响应变量的线性关系可以通过交互项的显著性判断
    interaction <- aov(weight ~ dose * gesttime, data = litter)
    summary(interaction)#交互项不显著
    
    
    #协方差分析  
    #必须为 响应变量 ~ 协方差 + 控制因素
    ancova <- aov(weight ~ gesttime +dose, data = litter)
    summary(ancova)#gesttime显著,说明对响应变量产生了影响
    #查看调整后的均值
    install.packages('effects')
    library(effects)
    effect('dose',ancova)
    
    #两两比较
    K <- rbind(contrMat(table(litter$dose), "Tukey")) 
                #根据‘Tukey’自动生成一个矩阵
    

    #Tukey为两两比较,有6种情况;Dunnet则是实验组与对照组比较,为3种(3行)

              
    mc1 <- glht(ancova, linfct = mcp(dose = K),#根据K矩阵计算,K可自定义
                alternative = "less")
    summary(mc1)
    
    #自定义两两比较
    K <- rbind('5-500'=c(0,1,0,-1),
               '0-500'=c(1,0,0,-1))
    
    mc2 <- glht(ancova, linfct = mcp(dose = K), alternative = "less")
    summary(mc2)
    
    ###根据不同方法调整P
    summary(mc1, test = univariate())
    summary(mc2, test = adjusted("bonferroni"))
    
    
    ##协变量对响应变量的线性关系可以通过交互项的显著性判断---之可视化判断####
    data1 <- subset(litter,dose==0)
    data2 <- subset(litter,dose==5)
    data3 <- subset(litter,dose==50)
    data4 <- subset(litter,dose==500)
    par(mfrow=c(2,2))
    plot(data1$gesttime,data1$weight,main = dose0)
    abline(lm(weight~gesttime,data=data1))
    
    plot(data2$gesttime,data2$weight,main = dose5)
    abline(lm(weight~gesttime,data=data2))
    
    plot(data3$gesttime,data3$weight,main = dose50)
    abline(lm(weight~gesttime,data=data3))
    
    plot(data4$gesttime,data4$weight,main = dose500)
    abline(lm(weight~gesttime,data=data4))

    多元方差

    #日照时间对农作物的 糖分、维生素、脂肪含量
    #课外阅读量对 学生语文、数学、历史成绩的影响
    #小学生每日VC摄入与其身高、胸围、BMI的影响
    
    
    #多元方差分析应用条件和要求:
    #多元方差分析指的是响应变量不止一个,且存在高度线性相关
    #对某个概念需要多个维度去进行描述,否则难以反应对象真实情况
    #要求自变量为分类变量
    
    rm(list = ls())
    
    library(MASS)
    data("UScereal")
    head(UScereal)
    #货架对各层营养含量的影响
    y <- with(UScereal, cbind(calories,fat,sugars))
    mav <- manova(y ~ shelf, data=UScereal)
    #两控制变量,3响应变量
    #mav <- manova(y ~ shelf + control +shelf*control, data=UScereal)
    
    summary(mav)
    ##更多前提假设参考
    #https://www.sohu.com/a/213290294_489312

    重复测量

    Valar morghulis
  • 相关阅读:
    4 Apr 18 软件开发目录 logging模块的使用 序列化(Json, Pickle) os模块
    3 Apr 18 内置函数 列表生成式与生成器表达式 模块的使用之import 模块的使用之from…import…
    2 Apr 18 三元表达式 函数递归 匿名函数 内置函数
    30 Mar 18 迭代器 生成器 面向过程的编程
    29 Mar 18 函数 有参、无参装饰器
    28 Mar 18 函数
    27 Mar 18 函数的参数
    26 Mar 18 函数介绍
    23 Mar 18 文件处理
    22 Mar 18 补充数据类型+字符编码+文件处理
  • 原文地址:https://www.cnblogs.com/super-yb/p/12733308.html
Copyright © 2011-2022 走看看