データを男女などで二つ以上の群にわけて、比較したいときは色々とあります。平均や分散の違いをt検定やF検定などで比較するのが教科書的かつ堅実な方法ですが、統計学に必ずしも詳しくない人々に説明する場合は、ビジュアライゼーションも考えていかないといけません。
この目的でよく使われるのは信頼区間によるエラーバーがついたグラフですが、実は往々にして誤解されています*1。2つの群のエラー・バーは重なっているのに、t検定で統計的に有意な差がでるときがあり、群を複数に分割したときの多重比較のための補正も行われません。比較対象が群ではなくて0などある一点であれば、信頼区間の厳密な説明*2とはちょっと違う使い方ではあるものの、t検定と合致した結果が得られるわけですが。
こういうわけで、群と群をグラフで比較するのには、信頼区間は使わない方が望ましいです。ガブリエル比較区間（Gabriel (1978)）と言う手法があるので、こちらを使うことをお勧めします。多重比較のための補正も行ってくれます。rgabrielと言うパッケージがあるので、使ってみましょう。
# 信頼区間のエラーバーが重なり、ガブリエル比較区間ではそうでないデータセット df1 <- data.frame(category = c("C","C","C","C","C","C","T-1","T-1","T-1","T-1","T-1","T-1"), value = c(0.117,-0.718,0.7199,-0.2422,-0.3391,-0.2223,1.0094,1.4306,0.099,1.2839,0.403,0.5616)) # 有意水準（片側） alpha <- 0.025 # ライブラリを呼び出す．インストールしていない場合は install.packages("rgabriel")を library(rgabriel) # 信頼区間とガブリエル比較区間を計算 df2 <- with(df1, { df2 <- data.frame( category=levels(as.factor(category)), count = as.numeric(tapply(value, category, length)), mean = as.numeric(tapply(value, category, mean)), sd = as.numeric(tapply(value, category, sd)) ) df2$se <- df2$sd/sqrt(df2$count) df2$t_upper <- df2$mean + df2$se*qt(1 - alpha, df=df2$count-1) df2$t_lower <- df2$mean - df2$se*qt(1 - alpha, df=df2$count-1) df2$g_upper <- df2$mean + rgabriel(df1$value, df1$category, alpha) df2$g_lower <- df2$mean - rgabriel(df1$value, df1$category, alpha) df2 })
データフレームdf2の中のt_upperとt_lowerが信頼区間、g_upperとg_lowerがガブリエル比較区間になります。信頼区間では群Cのt_upperは群T-1のt_lowerよりも大きいですが、g_upperはg_lowerよりも小さくなっています。グラフを描くと、以下のようにガブリエル比較区間の方が縮まっているのが分かります。
t検定をかけて、どちらが正しいかを見てみましょう。
t.test(df1$value[df1$category=="C"], df1$value[df1$category=="T-1"])
両側検定で5%有意で帰無仮説が棄却され、群Cと群T-1に差があることが分かります。なお、多重比較のための補正は手法によって結果が変わるので、検定結果と整合的になるとは限らないことには注意してください。
グラフのソースコード
gabriel.plotがあるので、plotすればグラフになるのですがビジュアル的にイマイチだったので、以下のように描画しました。
windowsFonts(Meiryo = windowsFont("Meiryo UI")) par(family="Meiryo", mar=c(2, 1, 1, 1), mgp=c(2, 1, 0)) unit <- 0.1 margin <- 0.05 # 信頼区間のエラーバーとガブリエル比較区間のエラーバーの隙間の幅 ylim <- c(unit*floor(min(df2$t_lower, df2$g_lower)/unit), unit*ceiling(max(df2$t_upper, df2$g_upper)/unit)) xlim <- c(c(0.5, length(df2$mean) + 0.5)) plot((1:length(df2$mean)) - margin/2, df2$mean, type="p", ylim=ylim, xlim=xlim, axes=FALSE, xlab="", ylab="", pch=19, col="red") arrows(1:length(df2$mean) - margin/2, df2$g_lower, 1:length(df2$mean) - margin/2, df2$g_upper, code=0, col="red", lwd=2) arrows(1:length(df2$mean) + margin, df2$t_lower, 1:length(df2$mean) + margin, df2$t_upper, code=0, lwd=2, lty=2) points(1:length(df2$mean) + margin, df2$mean, pch=19) abline(h=df2$t_upper[1], lty=2) abline(h=df2$g_upper[1], col="red", lty=3) abline(h=df2$t_lower[1], lty=2) abline(h=df2$g_lower[1], col="red", lty=3) axis(1, at=c(xlim[1], 1:length(df2$mean), xlim[2]), labels=c("",as.character(df2$category),""), pos=ylim[1], las=2) legend("topleft", c(sprintf("%.0f%%信頼区間", (1-2*alpha)*100), "ガブリエル比較区間"), col=c("black", "red"), lty=c(2, 1), lwd=c(2, 2), bg='white', box.col = 'black', bty="n")
メタアナリシス
rgabriel関数のソースコードを見ていると、群を表すファクターごとの標本標準偏差のベクトル（s）と、ファクターごとのサブサンプル・サイズのベクトル（n)と、有意水準（a）があればガブリエル比較区間は計算できるので、複数の分析のデータをつなぐことも簡単にできると思います。
…以上を書いて何年間も忘れていたのですが、rgabriel関数のコードを切り出して、標準誤差と自由度と有意水準から計算する関数を用意しました。
# se: 標準誤差のベクトル # df: 自由度のベクトル # a: 有意水準（片側） get_gci_by_se <- function(se, df, a=0.025){ if(length(se)!=length(df)){ return(NULL) } k <- length(se) # 属性の数 dfstar <- sum(df) # 全体の自由度 # the following code is written by Yihui according to Stoline's 1979 paper(equation 2.2). psmm_x = function(x, c, r, nu) { snu = sqrt(nu) sx = snu * x # for the scaled Chi distribution lgx = log(snu) - lgamma(nu/2) + (1 - nu/2) * log(2) + (nu - 1) * log(sx) + (-sx^2/2) exp(r * log(2 * pnorm(c * x) - 1) + lgx) } psmm = function(x, r, nu) { res = integrate(psmm_x, 0, Inf, c = x, r = r, nu = nu) res$value } qsmm = function(q, r, nu) { r = (r * (r - 1)/2) if (!is.finite(nu)) return(qnorm(1 - .5 *(1 - q^(1/r)))) res = uniroot(function(c, r, nu, q) { psmm(c, r = r, nu = nu) - q }, c(0, 100), r = r, nu = nu, q = q) res$root } SR <- qsmm(1-a, k, dfstar) # 元のソースのsをs = se/sqrt(n)に置き換えて整理 (vstar <- SR*se/sqrt(2)) }
