濃淡のある信用区間のプロット

WinBUGSやJAGSを用いたベイズ推定を使うと、各パラメーターに対する大量のMCMCサンプルが得られます。このMCMCサンプルを使って95%信用区間を計算し、それをもとに議論しますが、こうした任意の信用区間ではなく、事後確率がどのような分布をしているのかのほうを見せたい場合もあります。そんな時に便利なパッケージdenstripを紹介します。

濃淡のあるパラメーター推定結果の図示

ここでは、行列XにパラメーターのMCMCサンプルの入っているとしてみましょう。各列が各パラメーター、各行がMCMCサンプルに対応します。

X <- matrix(NA,nrow=1000,ncol=3)
for(i in 1:ncol(X)){
  X[,i] <- rnorm(1000,rnorm(1,2,2),runif(1,0.1,1))  
}
colnames(X) <- c("p1","p2","p3")
> head(X)
           p1        p2       p3
[1,] 2.527981 0.9171227 4.432790
[2,] 3.007233 2.9618068 4.689230
[3,] 1.891604 2.8134809 4.616570
[4,] 2.741595 0.8363671 4.741110
[5,] 2.602177 3.3400819 4.502327
[6,] 1.603801 2.5635456 4.356336

このデータを使って、事後分布と比例するグラデーション付きのパラメーター推定値を図示してみます。手順としては、まず空のプロットを呼び出し、そこにdenstrip関数をつかってストライプを図示していくことになります。

library(denstrip) #load package
plot(0,type="n",xlim=range(X),ylim=c(1,3)) # make plot region
for(i in 1:3) denstrip(X[,i],at=i,mticks=median(X[,i]),mcol="white")
# at = i is y-axis values density stribe to be plotted
# mticks is the x-position for median or similar
# mcol is the tick color to be used

出来上がりはこんな感じです。濃淡は、特に指定しなければdensity関数から求められる確率密度カーネルを基に図示されます。

また、以下のようにすれば見栄えも調節できます。

plot(0,type="n",xlim=range(X),ylim=c(1,3), ann=F, axes=F)
for(i in 1:3){
  denstrip(X[,i],at=i,mticks=quantile(X[,i],c(0.025,0.5,0.975)),
            mcol=c("black","white","black"))
}#adding 95% credible intervals
box(bty="l")
axis(1); axis(2,at=c(1,2,3),labels=colnames(X),las=2)

濃淡のある線形回帰の信用区間

パラメーターの推定結果に加え、線形回帰の推定結果についても濃淡をつけて図示することが可能です。基本的な手順は一緒ですが、denstrip関数の代わりにdensregion関数を使います(densregion関数についての詳細はこちらを参照)。まず、説明変数xについて、線形モデルによる予測値のMCMCサンプルが入った行列を用意します。先ほどと同じように、列が予測値と対応しており、行がMCMCサンプルの繰り返しになります。

> head(MCMC[,1:3])
     ypred[1] ypred[2] ypred[3]…
[1,] 0.336092 0.333215 0.330362
[2,] 0.384028 0.379466 0.374958
[3,] 0.469216 0.467234 0.465260
[4,] 0.259841 0.256886 0.253964
[5,] 0.350182 0.347407 0.344654
[6,] 0.597118 0.593395 0.589696
…

このデータをもとに、下記のように整形していきます。

est <- apply(MCMC,2,median) # median estimate
x <- seq(min(P),max(P),length=100) # range of explanatory variable
y <- seq(0,1.6,length=512) # range of y-axis; probability density will be estimated in this region
z <- matrix(nrow=length(x),ncol=length(y)) # empty case for density values
for(i in 1:length(x)) z[i,] <- density(MCMC[,i],from=min(y),to=max(y))$y # insert density values

library(denstrip)
plot(0,type="n",ylim=range(y),xlim=range(x), ann=F, axes=F)
densregion(x,y,z, colmax="red", gamma=0.8,scale=0.75,nlevels=75)
lines(x,est,lwd=2,col=grey(0,0.65))
axis(1);axis(2,las=2);box(bty="l")

こんな感じになります。

SPONSOR LINK

Spread the love