「真の割合の95%信頼区間」をブートストラップ法で推定するRスクリプト - 渋谷駅前で働くデータサイエンティストのブログ

「カイゼンしたらコンバージョン率が○○%→△△%にup!」は分母を無視したら成り立たないかもしれない

上記ブログ記事で話題に取り上げた「真の割合の95%信頼区間」の推定ですが、大したコードじゃないので一つのRスクリプト（関数）にまとめてみました。ソースはGitHubに置いてあります。

やってることはものすごく簡単で、要はサンプルサイズが小さい状況で○○%（例えば課金率とか翌週定着率とか）という「割合」の実測値がある場合に、それが統計学的に見て実際にはどれくらいばらつく可能性があるか？をシミュレーションするというものです。

そこで、以下の関数CIComputeに、真の割合の信頼区間を求めたいデータのサンプルサイズnpar、想定される割合nratio、ブートストラップ繰り返し計算回数nbootを与えれば、勝手に真の割合・標準偏差・95%信頼区間上限＆下限を返してくれて、ついでにブートストラップ標本のヒストグラムをプロットしてくれるようにしてあります。

CICompute <- function(npar,nratio,nboot) {
# npar: Sample size
# nratio: A ratio assumed as the real ratio
# nboot: The number of iteration for bootstrap resampling

  Data<-c(rep(1,ceiling(npar*nratio)),rep(0,floor(npar*(1-nratio))))
  Data<-as.matrix(Data)
  
  ResultBoot      <- matrix(NA, nboot, npar)
  
  for(i in 1:nboot){
    BootID          <- sample(1:nrow(Data), nrow(Data), replace = T)
    BootSample      <- Data[BootID, ]
    ResultBoot[i, ] <- mean(BootSample)
  }
  
  (MeanBoot <- mean(ResultBoot[, 1])) # Estimated ratio
  (SDBoot   <- sd(ResultBoot[, 1])) # Estimated SD of the ratio
  hist(ResultBoot) # Plotting a histogram of the bootstrap samples
  
  ci_low<-MeanBoot - 1.96*SDBoot # Lower bound of 95% CI
  ci_up<-MeanBoot + 1.96*SDBoot # Upper bound of 95% CI
  return(list(MeanBoot=MeanBoot,SDBoot=SDBoot,ci_low=ci_low,ci_up=ci_up,ResultBoot=ResultBoot))
}

例えば、以下のような感じで95%信頼区間の推定ができます。

> res<-CICompute(50,0.5,2000)
> res$MeanBoot
[1] 0.5045
> res$SDBoot
[1] 0.07063034
> res$ci_low
[1] 0.3660645
> res$ci_up
[1] 0.6429355
> res$ResultBoot
（以下略）

f:id:TJO:20130523154321p:plain

サンプルサイズがn=50だと、結構ばらつきますね。95%信頼区間が36.3-63.8%となってしまうので、うっかりすると真の割合が50%でも実測値が37%ぐらいになってしまうこともある、というわけです。

> res<-CICompute(200,0.8,2000)
> res$MeanBoot
[1] 0.805093
> res$SDBoot
[1] 0.02758377
> res$ci_low
[1] 0.7510288
> res$ci_up
[1] 0.8591571
> res$ResultBoot
（以下略）