【R】2群間においてt検定を実行し，アスタリスク付きのグラフを表示する

新しいよりよいグラフの記事が出ています

こちらをご覧ください。こっちのほうがきれい！！！！！！！！！！！！！

旧文書

1年前に自分用に作った2群間においてt検定を実行して棒グラフを作り，よくあるバーとアスタリスク(＊)をつけるRプログラムを載せておく．

これから卒論・修論の時期なので何か役に立つかもしれない．

two_barplot <- function(x, xlab = "", ylab = "", col = NA) {
    #testdata
    control <- c(3822,4103,4175,3736)
    treatment <- c(1063,1082,1385,1487)
    xy <- cbind(control, treatment)

    #もしCSVで読み込むなら上のブロックをコメントアウト
    #csv <- read.csv(x)
    #xy <- as.matrix(csv)

    #色の指定がなければ以下の色に設定する
    if (is.na(col)) {
    col <- c("#666666", "#CCCCCC")
    }
  
    #変数の設定
    sample.labels <- names(xy)
    pvalues <- rep(NA, length(sample.labels))
    
    #y軸の上限設定
    xmax  <- max(xy)
    xsd   <- sd(xy)
    
    #平均値及び標準偏差の算出
    xm <- apply(xy, 2, mean)
    xs <- apply(xy, 2, sd)
    
    #Pvalue の算出
    pvalues <- t.test(xy[,1], xy[,2])$p.value
    
    #graph size and fonts settings
    par(mar =c(5,15,5,15))
    par(family = 'sans')
    
    #棒グラフ描写
    b <- barplot(xm, beside = TRUE, col = col, border = col, xlab = xlab, ylab = ylab, ylim = c(0, (xmax + 1.1 * xsd)))
    
    #エラーバー描写
    arrows(b, xm-xs, b, xm+xs, code=3, lwd = 1, angle = 90, length = 0.1)
    
    #棒グラフを結ぶ線の点
    
    maxy = max(xm + xs) * 1.1
    stepy = max(xm +xs) * 0.1

    for (i in 1:length(pvalues)) {
      xi <- b[, i]
      yi <- xm + stepy * 1.5
      maxyi <- max(yi) + stepy
      if (pvalues[i] < 0.05) {
        lines(c(xi[1], xi[1], xi[2], xi[2]), c(yi[1], maxyi, maxyi, yi[2]))
        if (pvalues[i] < 0.001) {
          text((xi[1] + xi[2]) / 2, maxyi + stepy / 4, "***")
        } else if (pvalues[i] < 0.01) {
          text((xi[1] + xi[2]) / 2, maxyi + stepy / 4, "**")
        } else if (pvalues[i] < 0.05) {
          text((xi[1] + xi[2]) / 2, maxyi + stepy / 4, "*")
        }
      }
    }
    

    pvalues
    
}

上記を実行して関数を読み込ませた上で