はじめよう実験計画

実験を早く終わらせるための技術

phc曲線|『瀕死の統計学を救え』を参考に

phc曲線とは

先日『瀕死の統計学を救え』という本を読んで「phc曲線」なるものを勉強したので簡単にまとめてみたいと思います。

phc曲線とは、あるデータが得られたとき、その母数muの最尤推定を行い、その母数muがある値X以上(あるいは以下)である確率をプロットしたものです。

phc曲線を作図するためには、まず、実験で得られたデータの母数の最尤推定を行う必要があります。図1は最尤推定で得られる母数の事後分布の一例です。この事後分布はマルコフ連鎖モンテカルロ法という手法を使って作成されます。

f:id:Sturgeon:20200523112213p:plain

図1. 母数muの事後分布

図1に示したように、mu>2である領域がphc、つまり事後分布の累積確率分布がphc曲線になります。このようにして、ある値Xに対してphc(mu>X)をプロットしたphc曲線が図2です。phc曲線により、例えば母数muがX=2以上である確率が0.2(20%)であるというようんなことが読み取れます。

f:id:Sturgeon:20200523102644p:plain

図2. phc曲線の例

正規分布を仮定したときのphc曲線

それでは、実際にphc曲線を作図する具体的方法について書きます。2つの実験データ間(treatとcontrol)に差があるかどうかを表すphc曲線をRを使って描いてみました。

#マルコフ連鎖モンテルロ法のパッケージ
library(MCMCpack)
#乱数のセット
set.seed(1)
#実験群のデータ
treat <- c(5.5,10.0,5.5,6.0,9.0,9.5,6.5,7.0,12.5,7.0,6.5,10.5,9.0,4.5,6.5,9.5,10.0,9.5,10.5,6.5,8.5,12.5,4.0,9.0)
#コントロール群のデータ
control <- c(5.5,11.5,7.0,7.5,10.5,11.0,8.0,8.5,14.0,8.5,10.0,8.0,12.0,10.5,6.0,8.0,11.0,11.5,11.0,12.0,8.0,10.0,14.0,5.5,10.5,9.0)

#尤度関数。theta:推定する4つの母数(実験群の平均・標準偏差およびコントロール群の平均・分散)
#分布には正規分布を仮定
llnormfun <- function(theta, x) { a <- x[1:24] b <- x[25:50] if(theta[2] < 0.0 || treat[4] < 0){ l <- -Inf } else{ l <- sum(log(dnorm(a, mean = theta[1], sd=theta[2])))+ sum(log(dnorm(b, mean = theta[3], sd=theta[4]))) } return(l) } #事後分布f(theta|x) post.diff <- MCMCmetrop1R(llnormfun, theta.init=c(6,3,6,3), x=c(treat,control), thin=1, mcmc=10000, burnin=500, tune=1.0, verbose=500, logfun=TRUE) #事後分布のプロット plot(post.diff)

f:id:Sturgeon:20200522154305p:plain

図3. 4つの母数(実験群の平均・分散、コントロール群の平均・分散)の事後分布のプロット

 

続いて、以下のコードによりphc曲線を描きます。

#MCMCで生成したデータについて
#treatとcontrolの平均値の差
diff_of_mean <- post.diff[,3]-post.diff[,1]

#phc curve
phc <- function(v, hypo="larger", xlim=c(0,4), title=NULL){
  n <- 50
  X <- numeric(n)
  phc <- numeric(n)
  a <- seq(xlim[1],xlim[2],length=n)
  for (i in (1:n)) {
    X[i] <- a[i]
    if (hypo == "larger") {
      phc[i] <- length(v[v>a[i]])/length(v)  
    }
    else if (hypo == "smaller") {
      phc[i] <- length(v[abs(v)<a[i]])/length(v)  
    }
    else if (hypo == "non-overlap") {
      phc[i] <- length(v[abs(v)>a[i]])/length(v)  
    }
  }
  plot(X, phc, type = "b", col="blue", pch=16)
  title(title)
  grid()
}
#phc曲線(mu>X)
phc(diff_of_mean, "larger", title = "phc(mu>X)")
#phc曲線(|mu|<X)
phc(diff_of_mean, "smaller", title = "phc(|mu|<X)") 

図4にphc曲線を示します。今、mu = diff_of_meanなので、図3左は実験群とコントロール群の平均値の絶対値がXより大きい確率を、右図は平均値の絶対値がXより小さい確率を表しています。

例えば、図4左図において、実験群とコントロール群の差muが1.0より大きい確率は約0.7(70%)と読み取ることが出来ます。右図においては、muが1.0未満である確率が30%ということになります。

f:id:Sturgeon:20200522160307p:plain

図4. phc曲線。平均値の差の絶対値がXより大きい確率(左)、あるいは小さい確率(右)

二項分布を仮定したときのphc曲線

次に、2択問題に対してphc曲線を扱ってみました。正答数の分布に二項分布を仮定したときの成功確率pに関するphc曲線が、正答数(yes)と誤答数(no)の組合せを(1,1), (5,5), (50, 50)とした場合にどうなるか考えてみました。

#2択問題→YES or Noの配列
yes <- numeric(1)+1
no <- numeric(1)
yes5 <- numeric(5)+1
no5 <- numeric(5)
yes50 <- numeric(50)+1
no50 <- numeric(50)
#x: (yes, no)=(5,5), x10: (yes, no)=(10,10), x100: (yes, no)=(50,50) x <- c(yes, no) x <- sample(x) x10 <- c(yes5, no5) x10 <- sample(x10) x100 <- c(yes50, no50) x100 <- sample(x100) llnormfun <- function(p, x) { if(p < 0.0 || p > 1.0){ return(-Inf) } else{ l <- sum(log(dbinom(x, size = 1, p))) return(l) } } plot.samp <- MCMCmetrop1R(llnormfun, theta.init=0.8, x=x, thin=1, mcmc=10000, burnin=500, tune=1.0, verbose=500, logfun=TRUE) plot.samp10 <- MCMCmetrop1R(llnormfun, theta.init=0.8, x=x10, thin=1, mcmc=10000, burnin=500, tune=1.0, verbose=500, logfun=TRUE) plot.samp100 <- MCMCmetrop1R(llnormfun, theta.init=0.8, x=x100, thin=1, mcmc=10000, burnin=500, tune=1.0, verbose=500, logfun=TRUE) n <- 100 X <- numeric(n) phc <- numeric(n) phc4 <- numeric(n) phc100 <- numeric(n) a <- seq(0,1,length=n) for (i in (1:n)) { X[i] <- a[i] phc[i] <- length(plot.samp[plot.samp>a[i]])/length(plot.samp) phc4[i] <- length(plot.samp[plot.samp10>a[i]])/length(plot.samp) phc100[i] <- length(plot.samp[plot.samp100>a[i]])/length(plot.samp) print(i) } plot(X, phc, ylim=c(0,1.0), type = "b", col="red") par(new=T) plot(X, phc4, ylim=c(0,1.0), type = "b", col="blue", ylab = "") par(new=T) plot(X, phc100, ylim=c(0,1.0), type = "b", col="black", ylab = "")

図5に、正答数(yes)と誤答数(no)の組合せを(1,1), (5,5), (50, 50)とした場合のそれぞれについて、phc曲線を示しました。

縦軸の値phcは、正答率がX以上である確率です。たとえば、(yes, no)=(1,1)の場合(赤線)なら、正答率が0.8以上である確率は10%ということになります。

phc曲線は回答数が多くなるにつれて、なだらかなカーブから急なカーブに移っていくことがわかります。これは、回答数が多いほど、正答率の推定が正確になることに対応しています。例えば、(yes, no)=(50,50)のphc曲線(黒線)を見ると、この曲線が0<phc<1の値をとるのは0.4<X<0.6ですね。これは「正答率が0.4以上である確率が約100%、正答率が0.6以上である確率が約0%」ということを意味します。言い換えると

「正答率は0.4~0.6であるということが、ほぼ100%正しい」と判断することが出来るのです。

f:id:Sturgeon:20200522164659p:plain

図5. 2択問題の正答率に関するphc曲線

まとめ

本記事で参考にした『瀕死の統計学を救え』では、有意差検定の代替としてphc曲線の使用が推奨されています。個人的にも、ほとんど実用上の差がない場合に「統計的有意」と判断する有意差検定に納得のいかない部分がありました。phc曲線では、ある値以上の差があるかどうかの確実さを、曲線のなだらかさで示すことができる(図5に示すようにデータが増えるにしたがって、より多くの確証を持って結論を出せる)ので、有意差検定に替わる良い手段だなと思いました。