赖江山的博客分享 http://blog.sciencenet.cn/u/laijiangshan 生态、统计与R语言

博文

方差分析多重比较出图

已有 12763 次阅读 2018-12-3 22:37 |系统分类:科研笔记| 方差分析多重比较

下面这些代码将展示由Francois Gillet编写的两个方差分析多重比较的函数 boxplert()和boxplerk(),来源Numerical Ecology with R (second Edition)

#普通方差分析多重比较
boxplert <-  function(X,
           Y,
           main = NULL,
           xlab = NULL,
           ylab = NULL,
           bcol = "bisque",
           p.adj = "none",
           cexy = 1,
           varwidth = TRUE,
           las = 1,
           paired = FALSE)
  {#Y=factor(Y,levels=unique(Y)) #如果运行次代码,表示组别是按照原来的顺序,而不是字母顺序排列
    aa <- levels(as.factor(Y))
    an <- as.character(c(1:length(aa)))
    tt1 <- matrix(nrow = length(aa), ncol = 6)    
    for (i in 1:length(aa))
    {
      temp <- X[Y == aa[i]]
      tt1[i, 1] <- mean(temp, na.rm = TRUE)
      tt1[i, 2] <- sd(temp, na.rm = TRUE) / sqrt(length(temp))
      tt1[i, 3] <- sd(temp, na.rm = TRUE)
      tt1[i, 4] <- min(temp, na.rm = TRUE)
      tt1[i, 5] <- max(temp, na.rm = TRUE)
      tt1[i, 6] <- length(temp)
    }
    
    tt1 <- as.data.frame(tt1)
    row.names(tt1) <- aa
    colnames(tt1) <- c("mean", "se", "sd", "min", "max", "n")
    
    boxplot(
      X ~ Y,
      main = main,
      xlab = xlab,
      ylab = ylab,
      las = las,
      col = bcol,
      cex.axis = cexy,
      cex.lab = cexy,
      varwidth = varwidth
    )    
    require(agricolae)
    Yn <- factor(Y, labels = an)
    sig <- "ns"
    model <- aov(X ~ Yn)    
    if (paired == TRUE & length(aa) == 2)
    {
      coms <- t.test(X ~ Yn, paired = TRUE)
      pp <- coms$p.value
    }    else
    {
      pp <- anova(model)$Pr[1]
    }    
    if (pp <= 0.1)
      sig <- "."
    if (pp <= 0.05)
      sig <- "*"
    if (pp <= 0.01)
      sig <- "**"
    if (pp <= 0.001)
      sig <- "***"
    
    mtext(
      sig,
      side = 3,
      line = 0.5,
      adj = 0,
      cex = 2,
      font = 1
    )    
    if (pp <= 0.05) {
      comp <- LSD.test(model,                       "Yn",
                       alpha = 0.05,
                       p.adj = p.adj,
                       group = TRUE)      # gror <- comp$groups[order(comp$groups$groups), ]
      # tt1$cld <- gror$M
      gror <- comp$groups[order(rownames(comp$groups)), ]
      tt1$group <- gror$groups
      mtext(
        tt1$group,
        side = 3,
        at = c(1:length(aa)),
        line = 0.5,
        cex = 1,
        font = 4
      )
    }
    list(comparison = tt1, p.value = pp)
    
  }##非参数组间差异的Kruskal-Wallis检验
  boxplerk <-  function(X,
           Y,
           main = NULL,
           xlab = NULL,
           ylab = NULL,
           bcol = "bisque",
           p.adj = "none",
           cexy = 1,
           varwidth = TRUE,
           las = 1,
           paired = FALSE)
  {#Y=factor(Y,levels=unique(Y)) #如果运行次代码,表示组别是按照原来的顺序,而不是字母顺序排列
    aa <- levels(as.factor(Y))
    an <- as.character(c(1:length(aa)))
    tt1 <- matrix(nrow = length(aa), ncol = 7)    
    for (i in 1:length(aa))
    {
      temp <- X[Y == aa[i]]
      tt1[i, 1] <- mean(temp, na.rm = TRUE)
      tt1[i, 2] <- sd(temp, na.rm = TRUE) / sqrt(length(temp))
      tt1[i, 3] <- sd(temp, na.rm = TRUE)
      tt1[i, 4] <- min(temp, na.rm = TRUE)
      tt1[i, 5] <- max(temp, na.rm = TRUE)
      tt1[i, 6] <- median(temp, na.rm = TRUE)
      tt1[i, 7] <- length(temp)
    }
    
    tt1 <- as.data.frame(tt1)
    row.names(tt1) <- aa
    colnames(tt1) <- c("mean", "se", "sd", "min", "max", "median", "n")
    
    boxplot(
      X ~ Y,
      main = main,
      xlab = xlab,
      ylab = ylab,
      las = las,
      col = bcol,
      cex.axis = cexy,
      cex.lab = cexy,
      varwidth = varwidth
    )    
    require(agricolae)
    Yn <- factor(Y, labels = an)
    comp <- kruskal(X, Yn, p.adj = p.adj)
    sig <- "ns"
    
    if (paired == TRUE & length(aa) == 2)
    {
      coms <- wilcox.test(X ~ Yn, paired = TRUE)
      pp <- coms$p.value
    }    else
    {
      pp <- comp$statistics$p.chisq
    }    
    if (pp <= 0.1)
      sig <- "."
    if (pp <= 0.05)
      sig <- "*"
    if (pp <= 0.01)
      sig <- "**"
    if (pp <= 0.001)
      sig <- "***"
    
    gror <- comp$groups[order(rownames(comp$groups)), ]
    tt1$rank <- gror$X
    tt1$group <- gror$groups
    mtext(
      sig,
      side = 3,
      line = 0.5,
      adj = 0,
      cex = 2,
      font = 1
    )   
     if (pp <= 0.1)
      mtext(
        tt1$group,
        side = 3,
        at = c(1:length(aa)),
        line = 0.5,
        cex = 1,
        font = 4
      )
    
    list(comparison = tt1, p.value = pp)
    
  }
  # # 普通方差分析的多重比较
  library(agricolae)
## Warning: package 'agricolae' was built under R version 3.4.4
data(sweetpotato)# # 检验方差分析假设
shapiro.test(resid(aov(sweetpotato$yield ~ sweetpotato$virus)))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(aov(sweetpotato$yield ~ sweetpotato$virus))
## W = 0.95878, p-value = 0.7663
 bartlett.test(sweetpotato$yield, sweetpotato$virus)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  sweetpotato$yield and sweetpotato$virus
## Bartlett's K-squared = 2.3886, df = 3, p-value = 0.4958
boxplert(
   sweetpotato$yield,
  sweetpotato$virus,
  ylab = "yield",
   xlab = "virus",
   bcol = "orange",
   p.adj = "holm"
 )

## $comparison
##        mean       se       sd  min  max n group
## cc 24.40000 2.084067 3.609709 21.7 28.5 3     b
## fc 12.86667 1.246774 2.159475 10.6 14.9 3     c
## ff 36.33333 4.233727 7.333030 28.0 41.8 3     a
## oo 36.90000 2.482606 4.300000 32.1 40.4 3     a
## 
## $p.value
## [1] 0.0007334172

非参数Kruskal-Wallis检验多重比较

# # 案例
library(agricolae)
 library(stats)
 data(InsectSprays)
 boxplerk(
   InsectSprays$count,
   InsectSprays$spray,
   ylab = "count",
   xlab = "spray",
   bcol = "bisque",
  p.adj = "holm",
   las = 1)

## $comparison
##        mean        se       sd min max median  n     rank group
## A 14.500000 1.3623732 4.719399   7  23   14.0 12 52.16667     a
## B 15.333333 1.2329647 4.271115   7  21   16.5 12 54.83333     a
## C  2.083333 0.5701984 1.975225   0   7    1.5 12 11.45833     c
## D  4.916667 0.7225621 2.503028   2  12    5.0 12 25.58333     b
## E  3.500000 0.5000000 1.732051   1   6    3.0 12 19.33333    bc
## F 16.666667 1.7936476 6.213378   9  26   15.0 12 55.62500     a
## 
## $p.value
## [1] 1.510845e-10

方差分析多重比较.Rmd

请大家这里下载R markdown的文件



https://blog.sciencenet.cn/blog-267448-1149684.html

上一篇:NMDS的是否有解释率?
下一篇:R语言在生态学研究当中的应用分析
收藏 IP: 121.195.114.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (3 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-24 17:25

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部