p値の誤用がかなり多かったため、最近は点推定とp値ではなく、クリーブランドのドットプロットに、エラーバーのようなヒゲをつけた図で信頼区間を表示することが多くなりました。ひとつ真似して描いてみて、その問題点を確認の上、サブサンプルとサブサンプルの比較に便利なガブリエル比較区間のプロットを行ないます。
1 データセット
ランダム化比較実験を模して、介入群と対照群からなるデータセットを作ります。
df01 <- data.frame(
group = 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))2 点推定
群の種類とその数を把握しておきます。
clst <- levels(as.factor(df01$group))
noc <- length(clst)群ごとの平均値も取得しておきましょう。
means <- with(df01, tapply(value, group, mean))3 信頼区間
よく見かける信頼区間の図を描いて見ましょう。
3.1 計算
信頼区間を計算して行列にします。
ci <- sapply(1:noc, function(i){
with(df01, t.test(value[group==clst[i]]))$conf.int
})
dimnames(ci) <- list(c("lower", "upper"), clst)以下のような行列になります。
| C | T-1 | |
|---|---|---|
| lower | -0.6266025 | 0.2462922 |
| upper | 0.3983691 | 1.3495411 |
3.2 プロット
dotchartを呼んで、arrowsでヒゲを追記して終わりです。矢印部分の角度を90度にすることで、矢印らしくなくしています。
ヒゲが入りきるようにxlimには信頼区間の範囲を指定しています。
dotchart(as.numeric(means), xlim = 1.2*range(ci),
main = "95%信頼区間", cex = 1.25, labels = colnames(ci),
pch = 21, pt.cex = 1.5,
bg = "white", color = "black", lcolor = "darkgray")
arrows(ci[1, ], 1:noc, ci[2, ], 1:noc,
length = 0.1, angle = 90, code = 3)
4 信頼区間の利用上の注意点
信頼区間はあくまで定数より大や小を帰無仮説としたt検定との対応があるだけで、信頼区間と信頼区間は比較できません。
上述のデータセットも信頼区間を見ると領域が重なっているので、介入群と対照群に違いが無いように思えて来ますが、t検定を行なうと介入群と対照群が同じ分布から生成されたという帰無仮説は有意水準5%で棄却されます。
t.test(value ~ group, data=df01)
Welch Two Sample t-test
data: value by group
t = -3.1137, df = 9.9463, p-value = 0.01106
alternative hypothesis: true difference in means between group C and group T-1 is not equal to 0
95 percent confidence interval:
-1.5651545 -0.2589121
sample estimates:
mean in group C mean in group T-1
-0.1141167 0.7979167 また、信頼区間は多重比較補正を行なわないので、幾つも並べると誤った運用になります。
5 ガブリエル比較区間
介入群と対照群の比較をしたいのであれば、信頼区間ではなくガブリエル比較区間を使いましょう。二標本であればウェルチのt検定と一致した結果になりますし、多重比較補正もかかります。図と解釈の対応がよくなるので、プレゼンテーションにも優れています。
ガブリエル比較区間の計算は、rgabrielパッケージを用います。install.packages("rgabriel")などを行なって、インストールしてください。rgabrielが準備できたら、信頼区間のときと似たような手順で、計算、表示ができます。
library(rgabriel)
r_rgabriel <- rgabriel(df01$value, df01$group, 0.025)
gci <- sapply(1:noc, function(i){
c(means[i] - r_rgabriel[i], means[i] + r_rgabriel[i], use.names = FALSE)
})
dimnames(gci) <- list(c("lower", "upper"), clst)
dotchart(as.numeric(means), xlim = 1.2*range(gci),
main = "95%ガブリエル比較区間", cex = 1.25, labels = colnames(gci),
pch = 21, pt.cex = 1.5,
bg = "white", color = "black", lcolor = "darkgray")
arrows(gci[1, ], 1:noc, gci[2, ], 1:noc,
length = 0.1, angle = 90, code = 3)
6 信頼区間と比較区間の比較
信頼区間とガブリエル比較区間を並べてみましょう。
同じ値になるわけですが、信頼区間の上限と下限、ガブリエル比較区間の上限と下限がそれぞれあるので、平均値のデータを行列で渡すことにします。
c <- replicate(2, means)
colnames(c) <- c("CI", "GCI")行列の列が大項目、行が小項目になります。
| CI | GCI | |
|---|---|---|
| C | -0.1141167 | -0.1141167 |
| T-1 | 0.7979167 | 0.7979167 |
この行列を引数に、dotchartを呼び出します。
dotchart(c, xlim = 1.2*range(ci, gci), ylim = c(-1, 7),
main = "95%信頼区間とガブリエル比較区間", cex = 1.25,
pch = 21, pt.cex = 1.5,
bg = "white", color = "black", lcolor = "darkgray")
arrows(gci[1, ], 1:noc, gci[2, ], 1:noc,
length = 0.1, angle = 90, code = 3)
arrows(ci[1, ], 1:noc + nrow(gci) + 2, ci[2, ], 1:noc + nrow(gci) + 2,
length = 0.1, angle = 90, code = 3)
大項目と大項目は2つ縦軸の行間が空いているのに注意してください。