数学いらずの医科統計学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))
ここで、
・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)
95%CIのエラーバーが重ならないので、統計学的に差がある、ということになる。
論文でよくエラーバーが付いているが、あれはなんだ?
標準偏差（SD）をくっつけていることがある。これはデータ観測点の散らばりなので、重なっている/いないが何か結論をもたらしてくれるわけではない。
標本平均の標準誤差（SEM）をくっつけていることがある。これは標準偏差を標本数で割るため、標準偏差より幅は小さくなる。このため、もしエラーバーが重なるなら、p値はかなり大きい。
重ならない場合は、結論はよくわからない。
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)
エラーバーの話について、このような記事を教えて頂きました。
研究者の多くはエラーバーの意味をろくに理解していない
