Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2011-12-10

MikuHatsune2011-12-10

2つの平均値の比較:対応のないt検定

数学いらずの医科統計学PART6 CHAPTER30をRでやってみる。

対応のないt検定については、ただひたすらRだけをやる特別講義でも扱っている。

 

Frazier EP, et al (2006)は、神経伝達物質であるノルアドレナリンが、どの程度膀胱括約筋を弛緩させるか測定した。

old   <- c(20.8, 2.8, 50, 33.3, 29.4, 38.9, 29.4, 52.6, 14.3)
young <- c(45.5, 55, 60.7, 61.5, 61.1, 65.5, 42.9, 37.5)
boxplot(list(old, young))

f:id:MikuHatsune:20111210164036p:image

ここで、

・2つの平均値の差に対する信頼区間

・2つの平均値が同一であるとする帰無仮説を検定するp値

を考え、対応のないt検定について学ぶ。

t検定は、裏で正規分布を考えているパラメトリック手法である。ノンパラメトリック手法はいつか触れる。(以前やったことのあるブートストラップ法とは違う手法たち)

 

Rに放り込むとこうなる。ここでは、平均値の差が等しいという帰無仮説を、両側検定で考えている。

t.test(old, young, altenative="two.sided", conf.level=0.95)

	Welch Two Sample t-test

data:  old and young 
t = -3.6242, df = 13.778, p-value = 0.002828
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -37.501081  -9.590586 
sample estimates:
mean of x mean of y 
 30.16667  53.71250 

p値は0.003である。帰無仮説は棄却される。

 

2つの平均値の差に対する信頼区間を考える。

datamean <- sapply(list(old, young), mean)
stdev <- sqrt(unlist(lapply(list(old, young), var)))
n     <- unlist(lapply(list(old, young), length))
ciw   <- qt(0.975, n) * stdev / sqrt(n)
d.error    <- ciw   # 誤差範囲のデータ
d.error.x <- (0.2 + 1) * 1:length(datamean) - 1/2   # エラーバーを表示する位置(x座標)
barplot(datamean, names.arg=d.names, ylim=c(0,70), main="mean and 95%CI error bar")
arrows(d.error.x, datamean, d.error.x, datamean + d.error, angle=90)
arrows(d.error.x, datamean, d.error.x, datamean - d.error, angle=90)

f:id:MikuHatsune:20111210164037p:image

95%CIのエラーバーが重ならないので、統計学的に差がある、ということになる。

 

論文でよくエラーバーが付いているが、あれはなんだ?

標準偏差(SD)をくっつけていることがある。これはデータ観測点の散らばりなので、重なっている/いないが何か結論をもたらしてくれるわけではない。

 

標本平均の標準誤差(SEM)をくっつけていることがある。これは標準偏差を標本数で割るため、標準偏差より幅は小さくなる。このため、もしエラーバーが重なるなら、p値はかなり大きい。

重ならない場合は、結論はよくわからない。

f:id:MikuHatsune:20111210164038p:image

par(mfrow=c(1, 3))
# 95%信頼区間
datamean <- sapply(list(old, young), mean)
stdev <- sqrt(unlist(lapply(list(old, young), var)))
n     <- unlist(lapply(list(old, young), length))
ciw   <- qt(0.975, n) * stdev / sqrt(n)
d.error    <- ciw   # 誤差範囲のデータ
d.error.x <- (0.2 + 1) * 1:length(datamean) - 1/2   # エラーバーを表示する位置(x座標)
barplot(datamean, names.arg=d.names, ylim=c(0,70), main="mean and 95%CI error bar")
arrows(d.error.x, datamean, d.error.x, datamean + d.error, angle=90)
arrows(d.error.x, datamean, d.error.x, datamean - d.error, angle=90)

# SD
d.names   <- c("old", "young")
d.error   <- datasd   # 誤差範囲のデータ
d.error.x <- (0.2 + 1) * 1:length(datamean) - 1/2   # エラーバーを表示する位置(x座標)
barplot(datamean, names.arg=d.names, ylim=c(0,70), main="mean and sd error bar")
arrows(d.error.x, datamean, d.error.x, datamean + d.error, angle=90)
arrows(d.error.x, datamean, d.error.x, datamean - d.error, angle=90)

# SEM
datasem   <- unlist(lapply(list(old, young), sd))/unlist(lapply(list(old, young), length))
d.names   <- c("old", "young")
d.error   <- datasem   # 誤差範囲のデータ
d.error.x <- (0.2 + 1) * 1:length(datamean) - 1/2   # エラーバーを表示する位置(x座標)
barplot(datamean, names.arg=d.names, ylim=c(0,70), main="mean and sem error bar")
arrows(d.error.x, datamean, d.error.x, datamean + d.error, angle=90)
arrows(d.error.x, datamean, d.error.x, datamean - d.error, angle=90)

 

エラーバーの話について、このような記事を教えて頂きました。

研究者の多くはエラーバーの意味をろくに理解していない

rbykrbyk 2011/12/10 18:40 はじめまして。

エラーバーが重ならない場合にどう捉えればいいかは、
以下の記事が参考になるかもしれません。

『研究者の多くはエラーバーの意味をろくに理解していない』
http://d.hatena.ne.jp/kamedo2/20110224/1298536747

MikuHatsuneMikuHatsune 2011/12/11 17:06 ありがとうございます。参考としてリンクを貼らせていただきました。

ryamada22ryamada22 2011/12/13 05:56 視覚的情報から、「有意かどうか」の判断をさせると失敗する、という話しは面白いですね。
エラーバーを見せると、「有意かどうか」の判断(p値として捉える)には失敗しがちだ、という話しですね。
エラーバーという定義とp値という定義の2つの定義があって、それを計算式でつなぐことが要求されていて、難題だから失敗するのでしょうか。
エラーバーを使わずに「推定したものをもっときちんと(分布で)見せる」と、「2群は違うか」という判断には、失敗しないのではないか、という気もします。


エラーバーは、群ごとの推定値のはらつきの具合を視覚的に示している。
バーでないもので示すとすれば、ベル型の分布で示すことになる。
ベル型で示すとすると、2群のベル型分布の重なりを「見せる」こととなる。
その重なっている2つのベル型分布を見て、どのくらい重なっていたら、「めったに重なっていないことにするか」という話しとも言えるのでしょうか。
「片方の分布が重なる方向にあっ」て、かつ、「もう片方の分布も重なる方向にある」ときに、「2群の推定値には差がない」ことになります。
あることがおきて、かつ、別のことがおきる、という確率は、確率x確率で考えるので、95%CIを2つ置いて、それの重なりを考えると、100-95 % の出来事よりも珍しい→有意水準が厳しすぎる。

これは、2群を比べるデータでオッズ比を出すときとは少し違います、オッズ比(relative riskの推定値とみなしている)を計算して、その信頼区間を見て、それが1(帰無仮説)を含むかどうか、というときには、推定されているベル型の分布(オッズ比について推定した分布)が1個なので、2群のベル型分布2個での話しとは違った話しになるということでしょう

通りすがり通りすがり 2011/12/13 06:08 棒グラフの上にエラーバーを描くのは,「ダイナマイト・プロット」とか呼ばれ,良くないグラフとされているようです。
上に示されているページのように,描くのがよしとされるようですよ。

ryamada22ryamada22 2011/12/13 06:37 「ダイナマイト・プロット」と言うのは初めて聞きました。医学・生物学系論文でよく見ますが。ここにそれに関する日本語記事を見つけました:"http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc034/07521.html"。
"Dynamite plot"で検索すれば"http://emdbolker.wikidot.com/blog:dynamite","http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/DynamitePlots"
この箱ひげページ"http://flowingdata.com/2008/02/15/how-to-read-and-use-a-box-and-whisker-plot/"のコメントにもdynamite plotに関するコメントがありました。

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証