Rで棒グラフにエラーバーをつける

何かしらデータを取ったらグラフ作ってエラーバーを加えて
統計的有意差があるかを確認をするのだが
データ量が増えるとエクセルだとめんどい

Rで棒グラフを書いてエラーバーを作れるようにした

irisのデータでグラフを書く
結果はこんな感じ

元のデータこれ

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
~省略~

Rのコンソールで
> iris
で確認できる
アヤメ属に含まれる3種の花弁とがく片の測定データで
setosa, versicolor, virginicaそれぞれ50点ずつデータがある

まずデータを読み込み、種の名前を取り出す
種名はsetosa, versicolor, virginicaの3種
次にデータの項目
Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species
を取り出してグラフを書く関数に渡す

Data.Iris <- iris
Species.List <- as.character(unique(Data.Iris$Species))
Col.Data <- colnames(Data.Iris)
MakingBarGraph(Data.Iris, Species.List, Col.Data[5], Species.List[1])

関数MakingBarGraph()は引数に
irisのデータ、
種名のリスト、
種名の入ったカラムの名前、
統計的な比較をする種名、
を使う

今回はsetosaと他の2種の測定結果とに
統計的有意差があればスターをつける

エラーバー、スターの書き方は以下のサイトを参考にした
棒グラフ | カテゴリーで分けられるデータの視覚化（棒グラフおよびエラーバー）
出力はpng形式でファイル名はbargraph.png

グラフを書く関数がこちら

MakingBarGraph <- function(data.all, sp.list, col.name, sp.name){
  tmp <- subset(data.all, data.all[[col.name]] == sp.name)
  tmp <- tmp[, colnames(data.all) != col.name]

  data.mean <- c()
  data.error <- c()

  for(i in sp.list){
    tmp <- subset(data.all, data.all[col.name] == i)
    tmp <- tmp[, colnames(data.all) != col.name]
    # 統計量を求める、列方向に計算
    # StatCalは別に定義した関数
    data.all.stat <- (apply(tmp, 2, StatCal))
    data.all.stat <- data.frame(t(data.all.stat))

    data.mean <- rbind(data.mean, data.all.stat$Mean)
    data.error <- rbind(data.error, data.all.stat$SEM)
  }
  colnames(data.mean) <- colnames(tmp)

  rownames(data.mean) <- sp.list
  colnames(data.error) <- paste(colnames(tmp), "error", sep = ".")

  print(data.mean)
  print(data.error)

  png.name = sprintf("bargraph.png")
  png(png.name, width = 800, height = 600)

  bar.mean <- barplot(data.mean,
    col=c("#9bbb5960", "#c0504d60","#4f81bd60"), # bar color
    #main = File.Name.First, # graph title
    ylim = c(0, max(data.mean + data.error) + 1),
    beside = T
  )

  legend("topleft", legend = sp.list,
    fill =c("#9bbb5960", "#c0504d60","#4f81bd60"),
    col = c("#9bbb5960", "#c0504d60","#4f81bd60"),
    border = c("#9bbb5960", "#c0504d60","#4f81bd60"),
    box.lwd = 0, box.lty = 0
  )

  arrows(bar.mean, data.mean - data.error, bar.mean, data.mean + data.error,
    code = 3, lwd = 1, angle = 90, length = 0.1
  )

  # 有意差の検定
  T12 <- c()
  T13 <- c()
  for(i in 1:4){
    t1 <- data.all[data.all[[col.name]] == sp.list[1],]
    t2 <- data.all[data.all[[col.name]] == sp.list[2],]
    t3 <- data.all[data.all[[col.name]] == sp.list[3],]

    # まずはF検定 var.test()
    if(var.test(t1[,i], t2[,i])$p.value < 0.05){
      F12 = F
    }else{
      F12 = T
    }
    if(var.test(t1[,i], t3[,i])$p.value < 0.05){
      F13 = F
    }else{
      F13 = T
    }

    # T検定, t.test()
    T12 <- append(T12, t.test(t1[,i], t2[,i], var.equal = F12)$p.value)
    T13 <- append(T13, t.test(t1[,i], t3[,i], var.equal = F13)$p.value)
  }

  # **, *をグラフにつける

  for (i in 1:length(T12)) {
    xi <- bar.mean[2, i] # x axis
    yi <- data.mean[2, i] + data.error[2, i] + 0.5
    if (T12[i] < 0.05) {
      if (T12[i] < 0.01) {
        text(xi, yi, "**")
      } else if (T12[i] < 0.05) {
        text(xi, yi, "*")
      }
    }
  }
  # **, *をグラフにつける
  for (i in 1:length(T13)) {
    xi <- bar.mean[3, i] # x axis
    yi <- data.mean[3, i] + data.error[3, i] + 0.5
    if (T13[i] < 0.05) {
      if (T13[i] < 0.01) {
        text(xi, yi, "**")
      } else if (T13[i] < 0.05) {
        text(xi, yi, "*")
      }
    }
  }
  dev.off()
}

基本統計量はsummary()で求められるが
SEMは出てこないので基本統計量を求める関数をつかう
この回答を参考にした
How can I summarizing data statistics using R - Stack Overflow

StatCal <- function(x, Na.Remove=TRUE){
  result <- c(Mean = mean(x, na.rm=Na.Remove),
              SD = sd(x, na.rm=Na.Remove),
              Median = median(x, na.rm=Na.Remove),
              Min = min(x, na.rm=Na.Remove),
              Max = max(x, na.rm=Na.Remove),
              N = length(x),
              SEM = sd(x, na.rm=Na.Remove)/sqrt(length(x)))
}

もし規格化するなら次の関数を使う
呼び出し方はこんな感じ
Data.Iris <- NormData(Data.Iris, Col.Data[5], Species.List[1])
この場合はsetosaの平均値で規格化する
19 Septembe 2020 追記
平均値で割っているだけなので規格化とは異なるかも
normalizationと呼んでいる場合もあるようだけど
Normalizations: dividing by mean - Cross Validated
How, When, and Why Should You Normalize / Standardize / Rescale… – Towards AI

追記ここまで

NormData <- function(data.all, col.name, sp.name){
  # 列名の一時保存
  tmp.col <- data.all[, colnames(data.all) == col.name]
  tmp.col.name <- colnames(data.all)
  # 基準になるデータセットの取得と平均値の計算
  tmp <- subset(data.all, data.all[[col.name]] == sp.name)
  data.all.stat <- apply(tmp[, colnames(tmp) != col.name], 2, StatCal)
  data.all.stat <- data.frame(t(data.all.stat))
  data.all.mean <- data.all.stat$Mean
  # 規格化
  data.all <- t(apply(data.all[, colnames(data.all) != col.name], 1, function(x){x/data.all.mean}))
  data.all <- data.frame(data.all)
  data.all <- cbind(data.all, tmp.col)
  colnames(data.all) <- tmp.col.name

  return(data.all)
}

結果はこうなる

ggplotとか使えばもっと楽できれいにできるのかもだけど