誤差範囲

棒グラフにエラーバーを付ける方法

棒グラフにエラーバーを付ける方法。棒グラフは barplot 関数を利用する。また、エラーバーは arrows 関数を利用する。利用するサンプルデータは、イネに異なる濃度のジベレリン(AT、S5、S4、S3、S2)を与えてから一定期間を過ぎたイネ苗の丈を測ったものである。S4 群と S2 群には欠損値が存在するため、na.rm = TRUE オプションを利用して、欠損値を無視する。

x <- read.csv("http://stat.biopapyrus.net/data/tseed.csv",header = TRUE)
head(x)
##     At   S5   S4   S3   S2
## 1 1.29 2.25 3.91 4.59 5.58
## 2 0.82 2.31 4.24 4.95 4.26
## 3 0.87 2.40 2.81 4.78 4.19
## 4 1.12 2.49 4.25 4.29 5.78
## 5 1.09 2.58 3.75 5.11 4.10
## 6 0.99 2.28 3.49 3.51 5.20

m <- apply(x, 2, mean, na.rm = TRUE)  # 平均
s <- apply(x, 2, sd, na.rm = TRUE)    # 標準偏差

# 棒グラフを描く
# b には棒グラフの頂点座標が保存される
b <- barplot(m, ylim = c(0,6), ylab = "Length(mm)", xlab = "Group", names.arg=names(x))

# エラーバーを描く
arrows(b, m - s, b, m + s, angle = 90, length = 0.1)
arrows(b, m + s, b, m - s, angle = 90, length = 0.1)
棒グラフのエラーバー

折れ線グラフにエラーバーをつける方法

折れ線グラフにエラーバーを付ける方法。サンプルデータは、ラットに 3 種類の餌(PF、C、G)を与え、7 日間の体重を測ったものである。それぞれの餌群には 5 個体が存在する。

# データの読み込み
x <- read.table("http://stat.biopapyrus.net/data/rat.food.txt", header = TRUE)
head(x)
##    type Day0 Day1 Day2 Day3 Day4 Day5 Day6 Day7
## 1    PF  200  203  195  192  187  185  184  178
## 2    PF  200  187  175  165  183  178  176  172
## 3    PF  208  197  209  199  198  190  186  185
## 4    PF  208  194  202  197  194  192  192  187
## 5    PF  212  200  210  199  188  192  188  185
## 6     C  192  211  223  232  243  248  258  266


xaxis <- 1:ncol(x[, -1])          # x 軸座標を設定
cols <- c("red", "blue", "green") # PF を赤色、C を青色、G を緑色

# 準備
plot(0, 0, type = "n", xlim = range(xaxis), ylim = range(x[, -1]),
     xlab = "Date", ylab = "Weight")

# 平均値と標準偏差を計算して、書き加える
type <- unique(x[, 1])             # 餌の種類を取得
for (i in 1:length(type)) {
  m <- apply(x[x[, 1] == type[i], -1], 2, mean)
  s <- apply(x[x[, 1] == type[i], -1], 2, sd)
  points(xaxis, m, col = cols[i])
  lines(xaxis, m, col = cols[i])
  arrows(xaxis, m + s, xaxis, m - s, angle = 90, length = 0.1, col = cols[i])
  arrows(xaxis, m - s, xaxis, m + s, angle = 90, length = 0.1, col = cols[i])  
}

# グラフに凡例を書き入れる
legend("topleft", legend = type, pch = 1, lty = 1, col = cols)
折れ線グラフにエラーバーを書き入れます

棒グラフに有意差を示すアスタリスクをつける方法

手順としては、棒グラフおよびエラーがを書きれてから、lines 関数を利用して対照群と処理群の棒グラフを結びつけて、text 関数でアスタリスクを書き入れる。

アスタリスクを書き入れる関数を作成する。

biobarplot <- function(x, xlab = "", ylab = "", col = NA) {
    sample.labels <- names(x)
    condition.labels <- colnames(x$[[1]])

    if (is.na(col)) {
        col <- c("#666666", "#CCCCCC")
    }
    
    # 平均、標準偏差、p-value を保存する変数を用意する
    dfm <- dfs <- matrix(0, ncol = length(condition.labels), nrow = length(sample.labels))
    colnames(dfm) <- colnames(dfs) <- condition.labels
    rownames(dfm) <- rownames(dfs) <- sample.labels
    pvalues <- rep(NA, length(sample.labels))
    
    # 平均、標準偏差、p-value を計算する
    for (i in seq(x)) {
        dfm[i, ] <- apply(x[[i]], 2, mean, na.rm = TRUE)
        dfs[i, ] <- apply(x[[i]], 2, sd, na.rm = TRUE)
        x1 <- x[[i]][, 1]
        x2 <- x[[i]][, 2]
        pvalues[i] <- t.test(x1[!is.na(x1)], x2[!is.na(x2)])$p.value
    }
    
    # データ整形
    dfm <- t(dfm)
    dfs <- t(dfs)
    
    # y 座標およびアスタリスク関連のプロットの位置を調整するための変数を用意する
    maxy <- max(dfm + dfs) * 1.1
    stepy <- max(dfm + dfs) * 0.1

    # 棒グラフ
    bb <- barplot(dfm, beside = TRUE, ylim = c(0, maxy + 2 * stepy), col = col, border = col, xlab = xlab, ylab = ylab)
    
    # エラーバー
    arrows(bb, dfm - dfs, bb, dfm + dfs, angle = 90, length = 0.25 / length(sample.labels))
    arrows(bb, dfm + dfs, bb, dfm - dfs, angle = 90, length = 0.25 / length(sample.labels))
    
    # 有意差を示すアスタリスクを書き入れる(p < 0.01 ならば **;p < 0.05 ならば * とする。)
    for (i in 1:length(pvalues)) {
        xi <- bb[, i]
        yi <- dfm[, i] + 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.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, "*")
            }
        }
    }
    
    # グラフ凡例
    legend("topleft", legend = condition.labels, fill = col, col = col, border = col, box.lwd = 0, box.lty = 0)
}

データを与えて棒グラフを描く。

cell.a <- data.frame(Control = c(1.1, 1.2, 0.9),
                     Treated = c(1.5, 1.6, 1.3))
cell.b <- data.frame(Control = c(1.0, 0.8, 0.9),
                     Treated = c(1.8, 1.9, NA))
cell.c <- data.frame(Control = c(1.4, 1.5, 1.2, NA),
                     Treated = c(2.1, 2.0, 1.8, 1.9))

data <- list(`Cell A` = cell.a, `Cell B` = cell.b, `Cell C` = cell.c)

biobarplot(data, xlab = "Cell type", ylab = "Value")