某アニメのパクリだが、厳密には姉妹誌であるNature Methodsの話である。
Nat Methods. 2015 Feb 26;12(3):179-85.にp値についてわかりやすい話とおそらくRでやったと思われるシミュレーションがあるので再現してみる。
p値について理解は深まるし、どちらかというとRで図を再現するほうが時間がかかる。
これぐらいだったらオレもNature Medthodsに投稿できたんじゃね?
図1
統計学的仮説検定を行うわけだが、「答えが分かっている問題」に対して、確実な検算ができると理解が深まる。Rによるシミュレーションはもっともよい例である。
統計学的には、というか医学的にもそうだが、「真の値」というものが存在する、と仮定する。例えば、人の身長はみんななのだけれども、ばらつきというものを統計学も医学も許容するから、たくさん集めれば平均が
で、その周囲に山となる。これをシミュレーション的に作る。
もうひとつ、仮説検定をするためには比較されるべき相手が必要である。これの平均はであり、差が
となるような、分散1の分布を作るとこんな感じになる。
今の時点ではの差があるということがわかっているので、統計を持ち出す必要性はないけれども、なぜ統計学が必要かというと、この山の形がわからないからである。
m1 <- 0 # 1つめの分布の平均 m2 <- 0.5 # 比較する分布の平均 sd1 <- sd2 <- 1 # どちらも分散は同じにしておく alpha <- 0.05 # 有意水準 col0 <- c("black", grey(0.7)) # fig 1 x <- seq(-4, 4, length=1000) y <- cbind(dnorm(x, mean=m1, sd=sd1), dnorm(x, mean=m2, sd=sd2)) matplot(x, y, type="l", col=col0, lty=1, lwd=3)
図2
真の分布が分かっているけれども、実験上はわかっていないという設定で、サンプルサイズ10のデータを得たとして、2標本t検定をしよう。
よく聞くと思われる説明だと、t検定は平均を比較している。いま、帰無仮説は、AとBの(母集団の)平均は等しい、ということを仮定している。いまのシミュレーションでを知っている我々としては、(そんなバカな…)と思うが、実際に10と10のデータを手にしている状態では、そんな母集団のことなぞ知るはずがない。
というわけで、「2つの母集団平均が一緒だとしたときの、今回の10のデータを生じる確率」がp値として計算される。図では、AとBの10サンプルずつの平均と、2標本t検定のp値を計算しているが、平均の差が1以上になってp=0.001となったシミュレーションや、Bの方が平均値がAより小さくなってしまったシミュレーションも偶然存在している。
プロット用のjitter関数はここからパクる。
# fig 2 n <- 10 # サンプルサイズ niter <- 4 # シミュレーション回数 g <- rep(1:2, each=n) # 群 # ふたつの分布からランダムサンプリングする tmp <- replicate(niter, mapply(function(a, b) rnorm(n, a, b), list(m1, m2), list(sd1, sd2)))
# jitter 関数をパクる myBoxplotJitter <- function(x, range, overlap=1, lambda=5) { index <- seq(1, length(x)) jitter <- rep(0, length(x)) tie <- rep(FALSE, length(x)) df <- data.frame(x, jitter, index, tie) df <- df[order(df$x),] for (i in 1:length(x)) { region <- c(seq(i - overlap, i), seq(i + 1, i + overlap)) region <- region[region>0 & region<=length(x)] df$jitter[i] <- abs( sqrt(sum((df$x[region] - df$x[i])^2)) / (length(region) - 1) ) df$tie[i] <- (length(table(df$x[region])) < length(region)) } attenuation <- exp(-lambda*df$jitter) / exp(-lambda*min(df$jitter)) attenuation <- ifelse(df$tie, 1, attenuation) df$jitter <- runif(length(x), min=-abs(range), max=abs(range)) * attenuation df <- df[order(df$index),] df$jitter }
yl <- c(-2, 3) layout(matrix(seq(niter+1), nr=1), widths=c(0.4, rep(1, niter))) r <- 0.3 # jitter のばらつき par(mar=c(5, 0, 4, 0)) plot(1, type="n", axes=FALSE, xlab="", ylab="", ylim=yl) axis(2, line=-3) title("Estimated\neffect size") par(mar=c(5, 2, 4, 2)) for(i in seq(niter)){ x <- c(tmp[,,i]) es <- round(diff(tapply(x, g, mean)), 3) p <- round(t.test(x ~ g)$p.value, 3) jitter <- unlist(by(x, g, function(x) myBoxplotJitter(x, range=0.1, lambda=2))) plot(g + jitter, x, col=rep(col0, each=n), ylim=yl, pch=16, axes=FALSE, xlab="", ylab="", cex=2) title(paste("Simulation", i, "\n", es, paste("(P =", paste(p, ")", sep="")))) axis(1, at=1:2, labels=LETTERS[1:2], tick=FALSE, cex.axis=2) m12 <- tapply(x, g, mean) for(j in seq(2)) segments(j-r, m12[j], j+r, col=col0[j], lwd=2) }
図3
サンプルサイズが大きくなると、信頼区間の推定幅が狭くなる。
サンプルサイズが大きくなると、それだけ真の形に近づくということなので、あやふやさがなくなるので信頼区間が狭くなる。信頼区間とは、ここでのシミュレーションでは分布の真の平均値が含まれる区間のことである。
# fig 3 xl <- c(0, 1.6) niter <- 1000 # シミュレーション回数 ss <- c(10, 30, 64, 100) # sample size par(mfrow=c(1, 4), mar=c(5, 3, 4, 2), xpd=TRUE) for(i in seq(ss)){ tmp <- replicate(niter, mapply(function(a, b) rnorm(ss[i], a, b), list(m1, m2), list(sd1, sd2))) h <- mapply(function(x) diff(x$conf.int), apply(tmp, 3, t.test)) hist(h, xlim=xl, xlab="Width of confidence interval", ylab="", main="", yaxt="n") p <- par()$usr po <- round(power.t.test(ss[i], delta=m2-m1, sd=sd1)$power, 2) t0 <- paste("Sample size", ss[i], "\nTheoretical power", po) if(i == 1){ text(p[1], mean(p[3:4]), t0, pos=4, cex=1.5) } else { text(0.7, mean(p[3:4]), t0, pos=4, cex=1.5) } }
図4
サンプルサイズはp値の分布に影響する。
今回のシミュレーションでは、という差があることが既に分かっている。この状況でシミュレーションをしてt検定をすれば、小さなp値が出やすいのは分かっている。
ここで、サンプルサイズを大きくしていけば、図3でも言ったように、真の形に近づいていくので、帰無仮説：を信じるのは難しくなる。
の差があるということが分かっているとき(effect size という)、小さなp値が出やすくなり、統計学的有意差が出やすくなる。この統計学的有意差が出やすくなる確率は計算することができて、検出力という。各ヒストグラムの赤字は、p値が0.05を下回ったシミュレーションの回数を示しており、数学的にきちんと計算した検出力とほぼ等しくなっていることがわかる。
論文では触れられていないが、2群に全く差がないとき、つまり、というか2群と言わず、1つの分布からサンプリングを繰り返して2つのデータグループにしたときに、p値の分布はどうなるだろうか。答えは一様分布なのだが、時間があれば確かめて欲しい。
# fig 4 xl <- c(0, 1.6) niter <- 1000 ss <- c(10, 30, 64, 100) # sample size # ランダムサンプリングする tmp <- mapply(function(y) replicate(niter, mapply(function(a, b) rnorm(y, a, b), list(m1, m2), list(sd1, sd2))), ss) # two sample t.test をする p <- mapply(function(i) mapply(function(j) t.test(tmp[[i]][,1,j], tmp[[i]][,2,j])$p.value, seq(niter)), seq(ss)) # ヒストグラムを計算しておく h <- apply(log(p, 10), 2, hist, nclass=20) xl <- c(-4, 0) yl <- c(0, max(unlist(mapply(function(x) x$counts, h)))) pcut <- c(10^(-4:-2), alpha) par(mfrow=c(1, 4), mar=c(5, 3, 4, 2)) for(i in seq(ss)){ idx <- cumsum(h[[i]]$counts) < sum(p[,i] < alpha) col3 <- replace(rep("white", length(h[[i]]$counts)), idx, col0[2]) n_sig <- mapply(function(x) colSums(p <= x), pcut) plot(h[[i]], xlim=xl, ylim=yl, xlab=expression(log[10]*italic(P)), ylab="", main="", yaxt="n", col=col3) abline(v=log(pcut, 10), col=col0[2]) pa <- par()$usr po <- round(power.t.test(ss[i], delta=m2-m1, sd=sd1)$power, 2) text(log(pcut, 10), pa[4], n_sig[i,], pos=3, xpd=TRUE, col=replace(rep(1, length(ss)), pcut==alpha, "red"), cex=1.5) text(log(1, 10), pa[4], niter, pos=3, xpd=TRUE, cex=1.5) t0 <- paste("Sample size", ss[i], "\nTheoretical power", po) legend("topleft", legend=t0, bty="n", bg="white", cex=1.5) }
図5
サンプルサイズが小さいときの、統計学的有意な結果を得たときは、効果量を過大評価している。
いま、サンプルサイズをいくつかいろいろ設定して、いままでと同様な2標本t検定を行った。このとき、p値と、AとBの2群の平均の差(効果量)を計算している。
右の散布図はp値と効果量であり、左のヒストグラムは、であり、かつ、p値が0.05を下回ったシミュレーションでの効果量の分布である。
例えば、サンプルサイズが10のとき、検出力は18%しかないが、p値が0.05を下回った結果を得た時のは、真の値である0.5より上回る傾向がある。というか、p値が0.05を下回ったときはすべて
である。
例えば、サンプルサイズが64のとき、検出力は80%であるが、p値が0.05を下回った結果を得た時のは、真の値である0.5より上回る場合が66%存在している。
これは先に言ったように、サンプルサイズが少ないと推定精度が落ち、たくさんの差がないと統計学的有意差が得られにくいからである。そのため、標本平均が母集団平均に比べてものすごい差があるときにこのような現象が見られ、winner's curse と呼ばれる。
カーネル密度推定による二次元分布等高線はこちらやこちらをぱくる。ヒストグラムを横向きに書くにはこちらをパクる。
# fig 5 # tmp, p は共通 library(MASS) # 各シミュレーションでのサンプル平均の差をとる m0 <- mapply(function(y) apply(apply(tmp[[y]], 3, colMeans), 2, diff), seq(tmp)) xl <- c(-4, 0) yl <- c(0, max(m0)) # 平均の差が正で、統計学的有意差のある効果量だけ分布をとる h1 <- mapply(function(i) hist(m0[idx & p[,i] < alpha, i], nclass=20), seq(ss), SIMPLIFY=FALSE) par(mfrow=c(4, 2)) for(i in seq(ss)){ idx <- m0[,i] > 0 x0 <- log(p[idx, i], 10) y0 <- m0[idx, i] par(mar=c(2, 4, 2, 1.8)) tmp_m <- m0[idx & p[,i] < alpha, i] sig <- round(sum(tmp_m > m2-m1)/length(tmp_m), 2) A <- h1[[i]] colidx <- A$breaks[1:(length(A$breaks) - 1)] >= m2-m1 col4 <- replace(rep("white", length(colidx)), colidx, col0[2]) plot(NULL, type="n", xlab="", ylab="", xlim=rev(c(0, max(A$counts))), ylim=yl, axes=FALSE) rect(0, A$breaks[1:(length(A$breaks) - 1)], A$counts, A$breaks[2:length(A$breaks)], col=col4) abline(h=m2-m1, col="red") text(par()$usr[1], m2-m1, sig, pos=2, xpd=TRUE, cex=1.5) po <- round(power.t.test(ss[i], delta=m2-m1, sd=sd1)$power, 2) t0 <- paste("Sample size", ss[i], "\nTheoretical power", po) legend("topleft", legend=t0, bty="n", bg="white", cex=1.5) # カーネル密度推定をする kd <- kde2d(x0, y0, c(bandwidth.nrd(x0), bandwidth.nrd(y0)), n=1000) par(mar=c(2, 1, 2, 2)) #image(kd, xlim=xl, ylim=yl, col=grey(seq(1, 0, length=100)), las=1) #points(x0, y0, pch=16, col="red", cex=0.6) plot(x0, y0, xlab="", ylab="", pch=16, col="red", cex=0.6, xlim=xl, ylim=yl, las=1, cex.axis=1.5) abline(v=log(alpha, 10), h=m2-m1, lty=3, col=col0[2], lwd=2) contour(kd, xlab="latitude",ylab="longitude", add=TRUE) }
図6
サンプルサイズによって信頼区間の幅が変わる。
もう何回言ったことやら。
# fig 6 ss <- c(10, 30, 100) # サンプルサイズ niter <- 4 # シミュレーション回数 # ランダムサンプリング tmp <- mapply(function(y) replicate(niter, mapply(function(a, b) rnorm(y, a, b), list(m1, m2), list(sd1, sd2))), ss) # 効果量の推定 est <- mapply(function(i) mapply(function(j) diff(t.test(tmp[[i]][,1,j], tmp[[i]][,2,j])$estimate), seq(niter)), seq(ss)) # 信頼区間 CI <- matrix(mapply(function(i) mapply(function(j) -t.test(tmp[[i]][,1,j], tmp[[i]][,2,j])$conf.int, seq(niter)), seq(ss)), nr=2) x <- unlist(tmp) g <- unlist(mapply(rep, seq(2*length(ss)*niter), rep(ss, each=niter*2))) jitter <- unlist(by(x, g, function(x) myBoxplotJitter(x, range=0.1, lambda=2))) x_ss <- seq(ss)*niter*2+0.5 par(mfrow=c(2, 1), mar=c(1, 4, 3, 2)) plot(g + jitter, x, col=col0[(g+1)%%2+1], ylim=yl, pch=16, xlab="", ylab="Sample values", cex=0.7, axes=FALSE) axis(2, las=1) abline(h=0, lty=3) abline(v=x_ss) text((x_ss+c(0, head(x_ss, -1)))/2, par()$usr[4], ss, xpd=TRUE, pos=3) text(par()$usr[1], par()$usr[4], "Sample size", xpd=TRUE, pos=3) yl <- range(unlist(CI)) x_est <- colMeans(matrix(seq(length(ss)*niter*2), nr=2)) par(mar=c(1, 4, 1, 2)) plot(g + jitter, x, type="n", axes=FALSE, ylim=yl, xlab="", ylab="Effect size and 95% CIs") axis(2, las=1) abline(h=0) abline(h=m2-m1, lty=3) points(x_est, c(est), pch=16) for(i in seq(x_est)) mapply(function(j) arrows(x_est[i], est[i], y1=CI[j,i], angle=90, length=0.05), seq(2))
