ベイズ統計モデリングのまとめ

ここまで約1ヶ月間、集中的にベイズ統計について学習した。理解できたこととできていないことなど、未だ混沌としているが、「確率と情報の科学：データ解析のための統計モデリング入門（岩波書店、久保拓弥著）を中心に学習したことについて、これまでの知識も交えて以下に整理しておく。
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第１章ベイズ定理
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
ベイズ統計学は、1700年代にイギリスの牧師トーマス・ベイズにより発見されたという「ベイズの定理」に基づく。「ベイズの定理」とは、
事象Aが起こる確率をP(A)として、Aとは独立した事象Bが起こる確率をP(B)として場合、
事象Aが起きた後で、事象Bがおこる確率はP(B|A)であるが、P(B|A)にP(A)を掛け合わした確率P(B|A)*(A)（Aが起こる確率と、Aが起こったあとでBが起こる確率）は、その逆であるBが起こった後にAが起こる確率P(A|B)とBが起こる確率P(B)をかけ合わせたP(A|B)*P(B)と等しい。

P(B|A)*P(A) = P(A|B)*P(B)

この式を変形すればベイズの定理が得られる。

$P(B|A)={\frac {P(A|B)P(B)}{P(A)}}$

つまり、
P(B) = 事象 A が起きる前の、事象 B の確率（事前確率, prior probability）
P(B|A) = 事象 A が起きた後での、事象 B の確率（事後確率，条件付き確率, posterior probability，conditional probability）
とすると、ベイズの定理を使えば、事後確率 P(B|A) は上記に従って計算される。すなわち、事象Aに関するある結果（データ）が得られたとすると、それを反映し、事象Bが得られたあとのAが起こる確率P(A｜B)の乗算によって、事象 B の確率は事前確率から事後確率へと更新される。
ここでは、 P(A|B) は尤度(likelihood)と呼ばれる。

なお事象 B の確率の観点からは、P(A) は規格化定数としての意味しかないため、しばしば省略される。つまり事後確率は事前確率と尤度の積に比例する。

$P(B|A)\propto P(A|B)P(B)=P(A,B).$
————————————
ここで事象A＝原因として、事象Bを結果とすると
P(B) = 原因が起きる前の、結果の確率（事前確率, prior probability）
P(B|A) = 原因が起きた後での、結果の確率（事後確率，条件付き確率, posterior probability，conditional probability）

P(結果|原因) ∝ P(原因|結果)*P(結果)

原因から結果が起こる確率は、

事後確率 =　尤度(原因|結果)*事前確率

となる。具体的には、統計モデルに当てはめて推定事前分布を決定し、調査によって得られた尤度（結果が解っているときの原因の確率）を求め、両者の積として、原因に基づく結果の確率を算定する。
————————————
頭痛が起こる確率P(A)、インフルエンザに罹患する確率P(B)が解っているとする。
インフルエンザに罹患して頭痛が発生している人を調査したところ、P(A|B)が得られた。
頭痛を訴えた人の中でインフルエンザに罹患している確率P(B|A)は、

P(インフルエンザ|頭痛) = P(頭痛|インフルエンザ) * P(インフルエンザ) / P(頭痛)

もしくは

P(インフルエンザ|頭痛) ∝ P(頭痛|インフルエンザ) * P(インフルエンザ)

となる。
————————————
ベイズ統計は、ベイズの定理を確率分布に拡張し、

事後確率分布 =　尤度*事前確率分布
————————————

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第２章共役事前分布
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
統計学では、調査・実験によって得られたデータを元に、その未知なる母集団の分布なり、真の分布を推定して、未来に起こりうることに対する予測を得たいわけである。適切な予測をえるためには、無茶苦茶な分布ではなくて、なんらかのモデルにはまるような分布であってほしい。とくに数学的にハマる分布であれば、調査・実験結果を用いて、算術的に未来に対する予測を立てることが可能となる。

ベイズ定理により、事前分布を推定し、調査・実験によって得られた尤度をかけ合わせれば推定事後分布が得られる。性質がよく解っている事前分布を用いれば、「共役事前分布」を用いることで、複雑な計算を回避して簡便に推定事後分布が得られる。ここで、共役事前分布は、尤度をかけて事後分布を求めると、その関数形が同じになるような事前分布のことである。共役事前分布を使えば、事前分布と事後分布は同じ形になる。
性質がよく解っている事前分布と尤度、事後分布の関係は以下のようになっている。

共役事前分布 ? ? ?母数が規定する確率分布 ? ?事後分布
ベータ分布 ? ? ? ? ベルヌーイ分布　　　　 ? ?ベータ分布
ベータ分布 ? ? ? ? 二項分布　　　　　　　　ベータ分布
正規分布　　 ? ? 正規分布（σ2既知） ? ? ? ? ?正規分布
逆ガンマ分布 ? ? 正規分布（σ2未知） ? ? ? ? ? 逆ガンマ分布
ガンマ分布 ? ? ? ?ポアソン分布 ? ? ? ? ? ? ? ? ? ? ? ガンマ分布
ディリクレ分布 ? 多項分布 ? ? ? ? ? ? ? ? ? ? ? ? ? ?ディリクレ分布
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
このように共益事前分布を上手に利用して、調査・実験データから尤度を算定し、事前分布を推定して、ベイズの定理から事後分布を簡便に求めることができる。一方で、実際の統計データでは、常にこのような標準的な数学的分布にフィットする単純なモデルで描けるとは限らない。そこで、以下に示される一般線形モデルや、個体差を考慮した階層ベイズモデル（線形混合モデル）、さらに空間相関を考慮したCarモデルなどが有用となる。その場合、事後分布を推定する方法として、マルコフ・モンテカルロ法（MCMC法）が強力なツールとなる。
また、ベイズモデリングについて、BUGS言語を用いて策定し、BUGS言語に従ってベイズ解析を実行するWinBUGSをR言語ラッパー(R2WBwrapper.R)から操作する方法が利用できる。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第３章ベイズモデル化一般線形モデル　GLM:Generalized linear model
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

発生頻度がまれなる事象、例えば手術死亡率のようなものでれば、事前分布にポアソン分布を用いることが適当かもしれない。ポアソン分布の特徴は、平均値λが、分散σに等しい点である。

$P(X=k)={\frac {\lambda ^{k}e^{{-\lambda }}}{k!}}$

${\begin{aligned}\mathrm {E} [X]&=\lambda ,\\\mathrm {V} [X]&=\lambda .\end{aligned}}$

ここで、 $e$ はネイピア数 ( $e = 2.71828?$ )であり、 $k!$ は $k$ の階乗を表す。

また、 $λ$ は所与の区間内で発生する事象の期待発生回数に等しい。

$P (X = k)$ は、「所与の時間中に平均で $λ$ 回発生する事象がちょうど $k$ 回（ $k$ は非負の整数）発生する確率」に相当する。

発症数 yi は平均 λi のポアソン分布にしたがうとすると、

p(yi|λi)= λyi exp(−λi) / yi

i の平均 λi を以下のようにおいてみたら……

λi = exp(β1 + β2xi)

対数リンク関数

log(平均) = 線形予測子

線形予測子

?log(λi) = β1 + β2xi

β1は切片、β2は傾きとなる。

ポアソン分布は、平均(=分散）λで決定されるが、log(λi) = β1 + β2xiと規程したことで、

β1とβ2がポアソン分布を規定するパラメータとなる。

L(β1, β2) = Πp(yi|λi) = Πp(yi|β1,β2,xi)

X={xi}として、この尤度関数は、iで規程される離散確率分布なので、パラメーター{β1,β2}がある値を取っている時に

Yが得られる確率はp(Y|β1, β2) =L(β1, β2)

ベイズモデルでの事後分布は、（尤度）ｘ（事前分布）に比例するので、

p(β1, β2|Y) ∝?p(Y|β1,β2)p(β1)p(β2)

つまりp(β1, β2|Y)は、事後分布であり、データYが与えられたときのパラメーター{β1,β2}の同時確率である。

p(β1)とp(β2)は、それぞれ切片β1と傾きβ2の事前分布であり、これら事前分布を適切に指定すればベイズモデル化した一般線形モデルGLMとなる。

尤度は確率 (あるいは確率密度) の積であり，大量のかけ算となり、計算が大変なので、パラメーターの最尤推定では，対数尤度関数 (log likelihood function) を使う。対数尤度 log L(λ) の最大化は尤度 L(λ) の最大化になる。

ここで、パラメータβ1、β2の事前分布をどうするかである。事前に情報がない場合、ベイズ流では、無情報事前分布を選択する。ここでは、βの無情報事前分布として、平べったい正規分布N(0,100)を選ぶ。

> curve(dnorm(x, 0, 1.0E-4), -1, 1, lwd = 3)

1	> curve(dnorm(x, 0, 1.0E-4), -1, 1, lwd = 3)

dnorm(mean, tau) meanは平均、tauは分散の逆数, 標準偏差1/√tau
tau = 1/σ2=1.0E-4, σ^2 = 1000, σ= 100

サンプルデータをloadして、まずどんなデータか描画する。

> load("d.RData")
> d
          x  y
1  3.000000  5
2  3.210526  3
3  3.421053  6
4  3.631579  7
5  3.842105  7
6  4.052632  5
7  4.263158  9
8  4.473684  9
9  4.684211  7
10 4.894737 10
11 5.105263 12
12 5.315789  8
13 5.526316  7
14 5.736842  4
15 5.947368  4
16 6.157895 11
17 6.368421  9
18 6.578947  9
19 6.789474  8
20 7.000000  6

> load("d.RData")

> d

x y

1 3.000000 5

2 3.210526 3

3 3.421053 6

4 3.631579 7

5 3.842105 7

6 4.052632 5

7 4.263158 9

8 4.473684 9

9 4.684211 7

10 4.894737 10

11 5.105263 12

12 5.315789 8

13 5.526316 7

14 5.736842 4

15 5.947368 4

16 6.157895 11

17 6.368421 9

18 6.578947 9

19 6.789474 8

20 7.000000 6

xが体サイズで、yが種子数。観測データを図にしてみる。

> plot(d$x, d$y)

1	> plot(d$x, d$y)

サイズがxi=5の場合、種子数yiが10-12個がプロットされている。xi=5がピークのように見えるが、全体としてはばらついているように見える。

BUGSコードを用いて、このモデルを記述する：BUGSモデル：ポアソン分布

model
{
	Tau.noninformative <- 1.0E-4
	for (i in 1:N) {
		Y[i] ~ dpois(lambda[i])
		log(lambda[i]) <- beta1 + beta2 * (X[i] - Mean.X)
	}
	beta1 ~ dnorm(0, Tau.noninformative)
	beta2 ~ dnorm(0, Tau.noninformative)
}

model

{

Tau.noninformative <- 1.0E-4

for (i in 1:N) {

Y[i] ~ dpois(lambda[i])

log(lambda[i]) <- beta1 + beta2 * (X[i] - Mean.X)

}

beta1 ~ dnorm(0, Tau.noninformative)

beta2 ~ dnorm(0, Tau.noninformative)

}

次にBUGSの初期条件をRから入力していく。この際に、RからWinBUGSを操作するために「R2WBwrapper.R」を用いる。

> source("R2WBwrapper.R")
# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...
 要求されたパッケージ coda をロード中です 
 要求されたパッケージ boot をロード中です 
> load("d.RData")
> 
> clear.data.param()
> set.data("N", nrow(d))
> set.data("Y", d$y)
> set.data("X", d$x)
> set.data("Mean.X", mean(d$x))
> 
> set.param("beta1", 0)
> set.param("beta2", 0)

> ＃WinBUGSを呼び出し、サンプリング結果をpost.bugsに格納
> post.bugs <- call.bugs( 
+     file = "model.bug.txt",
+     n.iter = 1600, n.burnin = 100, n.thin = 3
+ )
> post.mcmc <- to.mcmc(post.bugs)

> source("R2WBwrapper.R")

# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...

要求されたパッケージ coda をロード中です

要求されたパッケージ boot をロード中です

> load("d.RData")

> clear.data.param()

> set.data("N", nrow(d))

> set.data("Y", d$y)

> set.data("X", d$x)

> set.data("Mean.X", mean(d$x))

> set.param("beta1", 0)

> set.param("beta2", 0)

> ＃WinBUGSを呼び出し、サンプリング結果をpost.bugsに格納

> post.bugs <- call.bugs(

+ file = "model.bug.txt",

+ n.iter = 1600, n.burnin = 100, n.thin = 3

+ )

> post.mcmc <- to.mcmc(post.bugs)

分析結果は、post.bugsに収まる。

R2WBwrapper.Rファイルで定義されているto.list()関数とto.mcmc()関数を用いて、post.bugsオブジェクトをそれぞれto.listクラスとmcmcクラスのオブジェクトに変換する。

> post.list <- to.list(post.bugs) 
> post.mcmc <- to.mcmc(post.bugs)
> s<-colnames(post.mcmc) %in% c("beta1", "beta2")
> s
[1]  TRUE  TRUE FALSE
> plot(post.list[,s,])

> post.list <- to.list(post.bugs)

> post.mcmc <- to.mcmc(post.bugs)

> s<-colnames(post.mcmc) %in% c("beta1", "beta2")

> s

[1] TRUE TRUE FALSE

> plot(post.list[,s,])

パラメータβ1とβ2の推定事後分布が以下のように得られた。

> plot(post.list[,"beta1",])

1	> plot(post.list[,"beta1",])

>plot(post.list[,"beta2",])

1	>plot(post.list[,"beta2",])

最後に、得られたパラメータβ1とβ2の推定事後分布を組み込んで、xiとyiの関係とλの事後分布のグラフを描く。
パラメータの事後分布からのサンプル

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> beta1 <- post.mcmc[, "beta1"]
> beta2 <- post.mcmc[, "beta2"]
> points(d$x, d$y)
> add.mean(median(beta1), median(beta2), lty = 1, lwd = 2)
> 
> par(ask = TRUE)
> plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
> for (i in 1:nrow(post.mcmc)) {
+     add.mean(beta1[i], beta2[i], lty = 1, col = "#00000004")
+ }
> points(d$x, d$y)
> add.mean(median(beta1), median(beta2), lty = 1, lwd = 2)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> beta1 <- post.mcmc[, "beta1"]

> beta2 <- post.mcmc[, "beta2"]

> points(d$x, d$y)

> add.mean(median(beta1), median(beta2), lty = 1, lwd = 2)

> par(ask = TRUE)

> plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")

> for (i in 1:nrow(post.mcmc)) {

+ add.mean(beta1[i], beta2[i], lty = 1, col = "#00000004")

+ }

> points(d$x, d$y)

> add.mean(median(beta1), median(beta2), lty = 1, lwd = 2)

パラメータβ1とβ2の関係も散布図に描く。

> plot(
+     as.matrix(post.mcmc)[,c("beta1", "beta2")],
+     lty = 1, col = "#00000030",
+     pch = 16, cex = 2,
+     xlab = "", ylab = ""
+ )

> plot(

+ as.matrix(post.mcmc)[,c("beta1", "beta2")],

+ lty = 1, col = "#00000030",

+ pch = 16, cex = 2,

+ xlab = "", ylab = ""

+ )

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第４章階層ベイズ化された一般線形混合モデル（個体差の有る場合）　GLMM:Generalized linear mixed model
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
次に、基本「有る無し」データのような二項分布に従うモデルで、個体差riを伴なう場合のモデル化を行う。
リンク関数と線形予測子を
logit(qi) = β + ri
とする。
切片β, 個体差riは、平均ゼロ標準偏差sの正規分布に従う。
p(ri|s)を階層事前分布と呼ぶ。「階層」とは、事前分布のパラメータにさらに事前分布が設定されていることをいう。また階層事前分布のパラメータsは超パラメータ（hyper parameter)、p(s)は超事前分布と呼ぶ。このような階層事前分布を使っているベイズ統計モデルが階層ベイズ統計モデルである。
ここでβには、無情報事前分布として平坦な正規分布N(0,100)を用いるとして、個体差riの階層事前分布には平均値ゼロで標準偏差がsの正規分布、この階層事前分布のばらつきsの事前分布は0

この階層ベイズモデルをBUGSで以下のように記述する。

model
{
for (i in 1:N) {
	Y[i] ~ dbin(q[i], 8) #二項分布
	logit(q[i]) <- beta + r[i] # 生存確率
}

解説：N個分ので観測データとの対応付け。観測された生存数Y[i]は、生存確率q[i]でサイズ8の二項分布(dbin(q[i]/ 8))に従う。種子の生存率q[i]は線形予測子beta + r[i]とロジットリンク関数logit(q[i])で与えられる。
以下は事前分布の定義

beta ~ dnorm(0, 1.0E-4) # 無情報事前分布>平均値ゼロで標準偏差10^2の正規分布
for (i in 1:N) {
  r[i] ~ dnorm(0, tau) # 階層事前分布>平均ゼロで標準偏差sの正規分布
}
tau <- 1 / (s * s)	 # tau tauは分散の逆数 
s ~ dunif(0, 1.0E+4) # 無情報事前分布＞0<s<10^4の一様分布
}

この階層ベイズモデルをBUGSで以下のように記述する。

model

{

for (i in 1:N) {

Y[i] ~ dbin(q[i], 8) #二項分布

logit(q[i]) <- beta + r[i] # 生存確率

}

解説：N個分ので観測データとの対応付け。観測された生存数Y[i]は、生存確率q[i]でサイズ8の二項分布(dbin(q[i]/ 8))に従う。種子の生存率q[i]は線形予測子beta + r[i]とロジットリンク関数logit(q[i])で与えられる。

以下は事前分布の定義

beta ~ dnorm(0, 1.0E-4) # 無情報事前分布>平均値ゼロで標準偏差10^2の正規分布

for (i in 1:N) {

r[i] ~ dnorm(0, tau) # 階層事前分布>平均ゼロで標準偏差sの正規分布

}

tau <- 1 / (s * s) # tau tauは分散の逆数

s ~ dunif(0, 1.0E+4) # 無情報事前分布＞0<s<10^4の一様分布

}

R2WBWrapperを用いて、BUGS言語で記載した統計モデルのパラメータに初期値を与えていく。

> source("R2WBwrapper.R")
# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...
> d <- read.csv("data7a.csv")
> 
> clear.data.param()
> set.data("N", nrow(d))
> set.data("Y", d$y)
> 
> set.param("beta", 0)
> set.param("r", rnorm(N, 0, 0.1))
> set.param("s", 1)
> set.param("q", NA)
NULL
> 
> post.bugs <- call.bugs(
+     file = "model.bug.txt",
+     n.iter = 10100, n.burnin = 100, n.thin = 10
+ )
> post.list <- to.list(post.bugs)
> post.mcmc <- to.mcmc(post.bugs)

> source("R2WBwrapper.R")

# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...

> d <- read.csv("data7a.csv")

> clear.data.param()

> set.data("N", nrow(d))

> set.data("Y", d$y)

> set.param("beta", 0)

> set.param("r", rnorm(N, 0, 0.1))

> set.param("s", 1)

> set.param("q", NA)

NULL

> post.bugs <- call.bugs(

+ file = "model.bug.txt",

+ n.iter = 10100, n.burnin = 100, n.thin = 10

+ )

> post.list <- to.list(post.bugs)

> post.mcmc <- to.mcmc(post.bugs)

得られたデータを元に、作図のRコードでグラフ表示する。

> d <- read.csv("data7a.csv")
> library(R2WinBUGS)
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> 
> q <- sum(d$y) / (8 * nrow(d))
> logistic <- function(z) 1 / (1 + exp(-z))
> n <- nrow(d)
> size <- 8
> q <- sum(d$y) / (size * n)
> 
> f.gaussian.binom <- function(alpha, x, size, fixed, sd)
+     dbinom(x, size, logistic(fixed + alpha)) * dnorm(alpha, 0, sd)
> 
> d.gaussian.binom <- function(v.x, size, fixed, sd) sapply(
+     v.x, function(x) integrate(
+         f = f.gaussian.binom,
+         lower = -sd * 10,
+         upper = sd * 10,
+         # for f.gaussian.binom
+         x = x,
+         size = size,
+         fixed = fixed,
+         sd = sd
+     )$value
+ )
> 
> plot.data <- function() plot(
+     0:size, summary(as.factor(d$y)),
+     ylim = range(c(0, dbinom(0:size, size, q) * n)),
+     xlab = "",
+     ylab = "",
+     pch = 19
+ )
> plot.polygon <- function(mm, p)
+ {
+     pp <- 1 - p
+     qp <- apply(mm, 1, quantile, probs = c(0.5 * pp, 1 - 0.5 * pp))
+     polygon(
+         c(0:size, size:0),
+         c(qp[1,], rev(qp[2,])),
+         border = NA,
+         col = "#00000020"
+     )
+ }
> 
> plot.lines <- function(mm)
+ {
+     apply(
+         mm, 2,
+         function(x) lines(0:size, x, col = "#00000001", lwd = 2)
+     )
+ }
> 
> beta <- post.mcmc[, "beta"]
> sigma <- post.mcmc[, "s"]
> if (!exists("my")) {
+     cat("# generating mp ")
+     mp <- sapply(
+         1:nrow(post.mcmc),
+         function(i) {
+             if (i %% 100 == 0) cat(".")
+             d.gaussian.binom(
+                 0:size, size,
+                 fixed = beta[i],
+                 sd = sigma[i]
+             )
+         }
+     )
+     cat(" done\n")
+     my <- apply(
+         mp, 2, function(prob) summary(
+             factor(
+                 sample(0:size, n, replace = TRUE, prob = prob),
+                 levels = 0:size
+             )
+         )
+     )
+ }
> plot.median <- function()
+ {
+     lines(
+         0:size, apply(mp * n, 1, median),
+         type = "b",
+         col = "black",
+         bg = "white",
+         pch = 21
+     )
+ }
> plot.data()
> q.env <- apply(my, 1, quantile, probs = c(0.025, 0.975))
> polygon(
+     c(0:size, size:0),
+     c(q.env[1,], rev(q.env[2,])), # 0.025 and 0.975
+     border = FALSE,
+     col = "#00000020"
+ )

> d <- read.csv("data7a.csv")

> library(R2WinBUGS)

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> q <- sum(d$y) / (8 * nrow(d))

> logistic <- function(z) 1 / (1 + exp(-z))

> n <- nrow(d)

> size <- 8

> q <- sum(d$y) / (size * n)

> f.gaussian.binom <- function(alpha, x, size, fixed, sd)

+ dbinom(x, size, logistic(fixed + alpha)) * dnorm(alpha, 0, sd)

> d.gaussian.binom <- function(v.x, size, fixed, sd) sapply(

+ v.x, function(x) integrate(

+ f = f.gaussian.binom,

+ lower = -sd * 10,

+ upper = sd * 10,

+ # for f.gaussian.binom

+ x = x,

+ size = size,

+ fixed = fixed,

+ sd = sd

+ )$value

+ )

> plot.data <- function() plot(

+ 0:size, summary(as.factor(d$y)),

+ ylim = range(c(0, dbinom(0:size, size, q) * n)),

+ xlab = "",

+ ylab = "",

+ pch = 19

+ )

> plot.polygon <- function(mm, p)

+ {

+ pp <- 1 - p

+ qp <- apply(mm, 1, quantile, probs = c(0.5 * pp, 1 - 0.5 * pp))

+ polygon(

+ c(0:size, size:0),

+ c(qp[1,], rev(qp[2,])),

+ border = NA,

+ col = "#00000020"

+ )

+ }

> plot.lines <- function(mm)

+ {

+ apply(

+ mm, 2,

+ function(x) lines(0:size, x, col = "#00000001", lwd = 2)

+ )

+ }

> beta <- post.mcmc[, "beta"]

> sigma <- post.mcmc[, "s"]

> if (!exists("my")) {

+ cat("# generating mp ")

+ mp <- sapply(

+ 1:nrow(post.mcmc),

+ function(i) {

+ if (i %% 100 == 0) cat(".")

+ d.gaussian.binom(

+ 0:size, size,

+ fixed = beta[i],

+ sd = sigma[i]

+ )

+ }

+ )

+ cat(" done\n")

+ my <- apply(

+ mp, 2, function(prob) summary(

+ factor(

+ sample(0:size, n, replace = TRUE, prob = prob),

+ levels = 0:size

+ )

+ }

> plot.median <- function()

+ {

+ lines(

+ 0:size, apply(mp * n, 1, median),

+ type = "b",

+ col = "black",

+ bg = "white",

+ pch = 21

+ )

+ }

> plot.data()

> q.env <- apply(my, 1, quantile, probs = c(0.025, 0.975))

> polygon(

+ c(0:size, size:0),

+ c(q.env[1,], rev(q.env[2,])), # 0.025 and 0.975

+ border = FALSE,

+ col = "#00000020"

+ )

> plot.median()

1	> plot.median()

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第５章　空間構造のある階層ベイズモデル
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
平均値λjを線形予測子と対数リンク関数を用いて、
logλj = β + ri
とする。
ここで空間の中で近傍のデータが影響するとする。
例えば、近傍とは空間的位置関係で隣の位置するものだけとし、その影響は均質に1/njだと規程する。
すると、ここでは近傍rj-1とrj+1の平均μiとし、μiが影響するとすれば、
rjの代わりにrj-μjを用いる正規分布

空間的な両端を除いて、
μi = (rj-1+rj+1)/k2
端の部分は、μ1=r2, μi=ri-1

p(rj|μj,s)　のばらつきのパラメータsはその場所でも同じと仮定する。このような事前分布は条件付き自己回帰(conditional auto regressive, CAR)モデルと呼ばれる。ここのでモデルのようにさまざまな制約をつけて簡単にしたものはintrinsic Gaussian CARモデルと分類される。

model
{
	for (j in 1:N.site) {
		Y[j] ~ dpois(mean[j])
		log(mean[j]) <- beta + r[j]
	}
	r[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau)
	beta ~ dnorm(0, 1.0E-4)
	tau <- 1 / (s * s)
	s ~ dunif(0, 1.0E+4)
}

model

{

for (j in 1:N.site) {

Y[j] ~ dpois(mean[j])

log(mean[j]) <- beta + r[j]

}

r[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau)

beta ~ dnorm(0, 1.0E-4)

tau <- 1 / (s * s)

s ~ dunif(0, 1.0E+4)

}

BUGSコードでの注意点は、
r[1:N.site]~car.normal(Adj[], Weights[], Num[], tau)

car.normal()関数は、intrinsic Gaussian CARモデルから[rj}のMCMCサンプルを発生させる関数である。Adj[]は隣接する場所の番号j、
Weights[]はAdj[]に対応する重み、
Num[]はそれぞれの区画に隣接する場所の数nj、
tauは分散の逆数、つまり1/s^2。

> load("Y.RData")
> set.data("N.site", length(Y))
> set.data("Y", Y)
> 
> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1)))
> set.data("Adj", c(2, Adj, N.site - 1))
> set.data("Weights", rep(1, 2 * N.site - 2))
> set.data("Num", c(1, rep(2, N.site - 2), 1))
> 
> set.param("beta", 0)
> set.param("r", rnorm(N.site, 0, 0.1))
> set.param("s", 1)
> 
> post.bugs <- call.bugs(
+     n.iter = 10100, n.burnin = 100, n.thin = 10
+ )
> post.list <- to.list(post.bugs)
> post.mcmc <- to.mcmc(post.bugs)

> load("Y.RData")

> set.data("N.site", length(Y))

> set.data("Y", Y)

> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1)))

> set.data("Adj", c(2, Adj, N.site - 1))

> set.data("Weights", rep(1, 2 * N.site - 2))

> set.data("Num", c(1, rep(2, N.site - 2), 1))

> set.param("beta", 0)

> set.param("r", rnorm(N.site, 0, 0.1))

> set.param("s", 1)

> post.bugs <- call.bugs(

+ n.iter = 10100, n.burnin = 100, n.thin = 10

+ )

> post.list <- to.list(post.bugs)

> post.mcmc <- to.mcmc(post.bugs)

> plot(Y, xlim=c(0,50), ylim=c(0,25))

1	> plot(Y, xlim=c(0,50), ylim=c(0,25))

グラフを描画させる。

> library(R2WinBUGS)
> if (!exists("Y")) load("Y.RData")
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> v <- 1:length(Y)
> y.max <- 27 
> 
> plot(
+     v, Y,
+     xlab = "location",
+     ylab = "y_i",
+     ylim = c(0, y.max)
+ )
> lines(m, lwd = 2, lty = 2)
> 
> mre <- sapply(
+     v, function(i) quantile(
+         post.mcmc[, sprintf("r[%i]", i)],
+         probs = c(0.5, 0.025, 0.975)
+     )
+ )
> b <- median(post.mcmc[,"beta"])
> polygon(
+     c(v, rev(v)),
+     exp(b + c(mre[2,], rev(mre[3,]))),
+     border = NA,
+     col = "#00000030"
+ )
> lines(exp(b + mre[1,]), lwd = 2)

> library(R2WinBUGS)

> if (!exists("Y")) load("Y.RData")

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> v <- 1:length(Y)

> y.max <- 27

> plot(

+ v, Y,

+ xlab = "location",

+ ylab = "y_i",

+ ylim = c(0, y.max)

+ )

> lines(m, lwd = 2, lty = 2)

> mre <- sapply(

+ v, function(i) quantile(

+ post.mcmc[, sprintf("r[%i]", i)],

+ probs = c(0.5, 0.025, 0.975)

+ )

> b <- median(post.mcmc[,"beta"])

> polygon(

+ c(v, rev(v)),

+ exp(b + c(mre[2,], rev(mre[3,]))),

+ border = NA,

+ col = "#00000030"

+ )

> lines(exp(b + mre[1,]), lwd = 2)

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
第6章　空間相関モデルと欠測のある観測データ
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
空間相関のモデルであるCarモデルは、GLMMモデルと比較して、観測データに欠損のあるデータ分析を適切にハンドリングできる。以下、欠損有るなし、GLMMかCarモデルかの4通りをまとめて推定して比較するBUGSモデルを作成する。

全部まとめて推定 (欠測あり・なし，CAR と GLMM)

model
{
	# j = 1: GLMM, no NA data　　GLMM 欠損データなし
	# j = 2: GLMM, NA data　　　 GLMM 欠損データあり
	# j = 3: CAR,  no NA data    CAR 欠損データなし
	# j = 4: CAR,  NA data　　　　CAR 欠損データあり
	Tau.noninformative <- 1.0E-4　　分散の逆数 Tau＝0.0001
	S.max <- 1.0E+4
	for (i in 1:N.site) {　　#i=1からi=50まで繰り返す
		for (j in 1:4) { #j=1からj=4まで繰り返す
			Y[i, j] ~ dpois(mean [i, j])　#Yはポアソン分布に従う。
			log(mean[i, j]) <- beta[j] + r[j, i]　#log(λ)=β + r[i]
		}
		for (j in 1:2) { 　#j=1からj=2繰り返す
			r[j, i] ~ dnorm(0.0, tau[j]) r[j,i]は正規分布N(0,分散の逆数)
		}
	}
	r[3, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[3])
        #r[3]について 
	r[4, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[4])
        #r[4]について
	for (j in 1:4) {　　j=1からj=4まで繰り返す
		beta[j] ~ dnorm(0, Tau.noninformative)　　#βは正規分布N(0, 100)
		tau[j] <- 1 / (s[j] * s[j])  #tauは、分散の逆数
		s[j] ~ dunif(0, S.max) #sには一様分布
	}
}

model

{

# j = 1: GLMM, no NA data　　GLMM 欠損データなし

# j = 2: GLMM, NA data　　　 GLMM 欠損データあり

# j = 3: CAR, no NA data CAR 欠損データなし

# j = 4: CAR, NA data　　　　CAR 欠損データあり

Tau.noninformative <- 1.0E-4　　分散の逆数 Tau＝0.0001

S.max <- 1.0E+4

for (i in 1:N.site) {　　#i=1からi=50まで繰り返す

for (j in 1:4) { #j=1からj=4まで繰り返す

Y[i, j] ~ dpois(mean [i, j])　#Yはポアソン分布に従う。

log(mean[i, j]) <- beta[j] + r[j, i]　#log(λ)=β + r[i]

}

for (j in 1:2) { 　#j=1からj=2繰り返す

r[j, i] ~ dnorm(0.0, tau[j]) r[j,i]は正規分布N(0,分散の逆数)

}

r[3, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[3])

#r[3]について

r[4, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[4])

#r[4]について

for (j in 1:4) {　　j=1からj=4まで繰り返す

beta[j] ~ dnorm(0, Tau.noninformative)　　#βは正規分布N(0, 100)

tau[j] <- 1 / (s[j] * s[j]) #tauは、分散の逆数

s[j] ~ dunif(0, S.max) #sには一様分布

}

post.bugsには、4種類(欠損なしGLMM、欠損ありGLMM、欠損なしCarモデル、欠損ありCarモデル)の分析結果βとr[i]が含まれている。
Yに対して欠損データYnoを作成するためのコードを以下に示す。
http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/set.data.R

> source("R2WBwrapper.R")
# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...
> clear.data.param()
> 
> load("Y.RData")
> no <- c(6, 9, 12, 13, 26:30)　#欠損データ番号をnoで指定する。
> Yno <- Y　　#データYをYnoへ
> Yno[no] <- NA　　#Ynoの[no]の位置をNA(欠損）とする。
> 
> set.data("N.site", length(Y))　#N.siteはYの数、ここでは50
> set.data("Y", matrix(c(Y, Yno, Y, Yno), N.site, 4)) 
　#欠損なし、欠損あり、欠損なし、欠損ありのデータYを行列Yに

> source("R2WBwrapper.R")

# reading "R2WBwrapper.R" (written by kubo@ees.hokudai.ac.jp)...

> clear.data.param()

> load("Y.RData")

> no <- c(6, 9, 12, 13, 26:30)　#欠損データ番号をnoで指定する。

> Yno <- Y　　#データYをYnoへ

> Yno[no] <- NA　　#Ynoの[no]の位置をNA(欠損）とする。

> set.data("N.site", length(Y))　#N.siteはYの数、ここでは50

> set.data("Y", matrix(c(Y, Yno, Y, Yno), N.site, 4))

　#欠損なし、欠損あり、欠損なし、欠損ありのデータYを行列Yに

欠損データあるなし、およびGLMMもしくはCarモデルの通りを上記BUGSモデルと、欠損データを作成するRを読み込んで、以下のように実行する。
http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/runbugs.R

source("set.data.R")　#set.data.Rを実行

set.param("beta", rep(0, 4))  #betaに{0,0,0,0}
set.param("r", matrix(rnorm(N.site * 4, 0, 0.1), 4, N.site))
set.param("s", rep(1, 4))#sに{1,1,1,1}

post.bugs <- call.bugs(n.iter = 10100, n.burnin = 100, n.thin = 10)
post.list <- to.list(post.bugs)
post.mcmc <- to.mcmc(post.bugs)

source("set.data.R")　#set.data.Rを実行

set.param("beta", rep(0, 4)) #betaに{0,0,0,0}

set.param("r", matrix(rnorm(N.site * 4, 0, 0.1), 4, N.site))

set.param("s", rep(1, 4))#sに{1,1,1,1}

post.bugs <- call.bugs(n.iter = 10100, n.burnin = 100, n.thin = 10)

post.list <- to.list(post.bugs)

post.mcmc <- to.mcmc(post.bugs)

ここでは、
post.bugs <- call.bugs(n.iter = 10100, n.burnin = 100, n.thin = 10) で指定されているように、10100回のサンプリング、burn-inの100回、10個ごとの間引きで、1ループで1000個のデータが生成される。デフォルトで3ループ繰り返されるので、3000データが生成されて、post.listには1000個づつｘ３ループ、post.mcmcには3000行のデータｘパラメータ数列となる。 post.mcmcのr[i]のデータからquantile()関数を用いて、mre[1]に中位点、mre[2]に2.5%位点、mre[3]に97.5%位点を求めている。今回の分析では、４つの異なる条件（欠損データあるなし、GLMMかCarモデルか）がまとめてpost.bugsに収められており、個体差のデータr[i]も二次元r[j,i]となっているために"r[%i]"を"r[j,%i]とする。また欠損なしのデータYに対して、欠損ありデータYnoを指定し、またパラメータbetaには、beta[1]?beta[4]を指定する。 1. 欠損なし、GLMMモデル

> library(R2WinBUGS)
> if (!exists("Y")) load("Y.RData")
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> v <- 1:length(Y)
> y.max <- 27 
> 
> plot(
+     v, Y,
+     xlab = "location",
+     ylab = "y_i",
+     ylim = c(0, y.max)
+ )
Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)
> 
> mre <- sapply(
+     v, function(i) quantile(
+         post.mcmc[, sprintf("r[1,%i]", i)],
+         probs = c(0.5, 0.025, 0.975)
+     )
+ )
> b <- median(post.mcmc[,"beta[1]"])
> polygon(
+     c(v, rev(v)),
+     exp(b + c(mre[2,], rev(mre[3,]))),
+     border = NA,
+     col = "#00000030"
+ )
> lines(exp(b + mre[1,]), lwd = 2)

> library(R2WinBUGS)

> if (!exists("Y")) load("Y.RData")

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> v <- 1:length(Y)

> y.max <- 27

> plot(

+ v, Y,

+ xlab = "location",

+ ylab = "y_i",

+ ylim = c(0, y.max)

+ )

Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)

> mre <- sapply(

+ v, function(i) quantile(

+ post.mcmc[, sprintf("r[1,%i]", i)],

+ probs = c(0.5, 0.025, 0.975)

+ )

> b <- median(post.mcmc[,"beta[1]"])

> polygon(

+ c(v, rev(v)),

+ exp(b + c(mre[2,], rev(mre[3,]))),

+ border = NA,

+ col = "#00000030"

+ )

> lines(exp(b + mre[1,]), lwd = 2)

2.欠損あり、GLMMモデル

> library(R2WinBUGS)
> if (!exists("Y")) load("Y.RData")
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> v <- 1:length(Y)
> y.max <- 27 
> 
> plot(
+     v, Yno,
+     xlab = "location",
+     ylab = "y_i",
+     ylim = c(0, y.max)
+ )
Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)
> 
> mre <- sapply(
+     v, function(i) quantile(
+         post.mcmc[, sprintf("r[2,%i]", i)],
+         probs = c(0.5, 0.025, 0.975)
+     )
+ )
> b <- median(post.mcmc[,"beta[2]"])
> polygon(
+     c(v, rev(v)),
+     exp(b + c(mre[2,], rev(mre[3,]))),
+     border = NA,
+     col = "#00000030"
+ )
> lines(exp(b + mre[1,]), lwd = 2)

> library(R2WinBUGS)

> if (!exists("Y")) load("Y.RData")

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> v <- 1:length(Y)

> y.max <- 27

> plot(

+ v, Yno,

+ xlab = "location",

+ ylab = "y_i",

+ ylim = c(0, y.max)

+ )

Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)

> mre <- sapply(

+ v, function(i) quantile(

+ post.mcmc[, sprintf("r[2,%i]", i)],

+ probs = c(0.5, 0.025, 0.975)

+ )

> b <- median(post.mcmc[,"beta[2]"])

> polygon(

+ c(v, rev(v)),

+ exp(b + c(mre[2,], rev(mre[3,]))),

+ border = NA,

+ col = "#00000030"

+ )

> lines(exp(b + mre[1,]), lwd = 2)

3. 欠損なし、Carモデル

> library(R2WinBUGS)
> if (!exists("Y")) load("Y.RData")
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> v <- 1:length(Y)
> y.max <- 27 
> 
> plot(
+     v, Y,
+     xlab = "location",
+     ylab = "y_i",
+     ylim = c(0, y.max)
+ )
Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)
> 
> mre <- sapply(
+     v, function(i) quantile(
+         post.mcmc[, sprintf("r[3,%i]", i)],
+         probs = c(0.5, 0.025, 0.975)
+     )
+ )
> b <- median(post.mcmc[,"beta[3]"])
> polygon(
+     c(v, rev(v)),
+     exp(b + c(mre[2,], rev(mre[3,]))),
+     border = NA,
+     col = "#00000030"
+ )
> lines(exp(b + mre[1,]), lwd = 2)

> library(R2WinBUGS)

> if (!exists("Y")) load("Y.RData")

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> v <- 1:length(Y)

> y.max <- 27

> plot(

+ v, Y,

+ xlab = "location",

+ ylab = "y_i",

+ ylim = c(0, y.max)

+ )

Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)

> mre <- sapply(

+ v, function(i) quantile(

+ post.mcmc[, sprintf("r[3,%i]", i)],

+ probs = c(0.5, 0.025, 0.975)

+ )

> b <- median(post.mcmc[,"beta[3]"])

> polygon(

+ c(v, rev(v)),

+ exp(b + c(mre[2,], rev(mre[3,]))),

+ border = NA,

+ col = "#00000030"

+ )

> lines(exp(b + mre[1,]), lwd = 2)

4. 欠損あり、Carモデル

> library(R2WinBUGS)
> if (!exists("Y")) load("Y.RData")
> if (!exists("post.bugs")) load("post.bugs.RData")
> #post.mcmc <- to.mcmc(post.bugs)
> post.mcmc <- as.mcmc(post.bugs$sims.matrix)
> v <- 1:length(Y)
> y.max <- 27 
> 
> plot(
+     v, Yno,
+     xlab = "location",
+     ylab = "y_i",
+     ylim = c(0, y.max)
+ )
Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)
> 
> mre <- sapply(
+     v, function(i) quantile(
+         post.mcmc[, sprintf("r[4,%i]", i)],
+         probs = c(0.5, 0.025, 0.975)
+     )
+ )
> b <- median(post.mcmc[,"beta[4]"])
> polygon(
+     c(v, rev(v)),
+     exp(b + c(mre[2,], rev(mre[4,]))),
+     border = NA,
+     col = "#00000030"
+ )
> lines(exp(b + mre[1,]), lwd = 2)

> library(R2WinBUGS)

> if (!exists("Y")) load("Y.RData")

> if (!exists("post.bugs")) load("post.bugs.RData")

> #post.mcmc <- to.mcmc(post.bugs)

> post.mcmc <- as.mcmc(post.bugs$sims.matrix)

> v <- 1:length(Y)

> y.max <- 27

> plot(

+ v, Yno,

+ xlab = "location",

+ ylab = "y_i",

+ ylim = c(0, y.max)

+ )

Hit <Return> to see next plot: lines(m, lwd = 2, lty = 2)

> mre <- sapply(

+ v, function(i) quantile(

+ post.mcmc[, sprintf("r[4,%i]", i)],

+ probs = c(0.5, 0.025, 0.975)

+ )

> b <- median(post.mcmc[,"beta[4]"])

> polygon(

+ c(v, rev(v)),

+ exp(b + c(mre[2,], rev(mre[4,]))),

+ border = NA,

+ col = "#00000030"

+ )

> lines(exp(b + mre[1,]), lwd = 2)

このように、Carモデルが欠損データを上手に取り扱えることが示された。

Science To Medicine

Just My Daily Study Note by ts.anesth.kpum

ベイズ統計モデリングのまとめ