準備

今回の実習に必要なパッケージを（必要ならまずインストールし、）呼び出す。

library("ggplot2")
library("plyr")
library("rgl")

多変量モデル：BDA3 (pp.74-78) 3.7 生物学的検定実験

データ

BDA3, p.74 のデータ (Racine et al. 1986) を使う。

Bio <- data.frame(x = c(-0.86, -0.3, -0.05, 0.73),
                  n = rep(5, 4),
                  y = c(0, 1, 3, 5))

推定するモデル

実験\(i = 1, \dots, k\) における薬物の投与量 \(x_i\)（単位はlog g/ml)、 動物の個体数 \(n_i\)、死亡した個体数\(y_i\) とする。

核実験内において、各個体が死ぬかどうか独立だと仮定すると、 \[y_i | \theta_i \sim \mathrm{Bin}(n_i, \theta_i)\] となる。ここで、\(\theta_i\) は死亡確率であり、これが推定の対象である。

死亡確率は薬物の投与量によって変化すると想定されるので、これをモデル化する。 ここでは、ロジスティック回帰モデルを使い、 \[\mathrm{logit}(\theta_i) = \alpha + \beta x_i\] あるいは \[\theta_i = \mathrm{logit}^{-1}(\alpha + \beta x_i) = \frac{1}{1 + \exp(-\alpha - \beta x_i)}\] とする。 したがって、このモデルでは、\(\theta_i\) をいう1つの変量を推定するのではなく、 \(\alpha\) と \(\beta\) という2変量を推定することになる。

尤度、事前分布、事後分布

このモデルの尤度は、 \[p(y_i | \alpha, \beta, n_i, x_i) = \theta_i^{y_i} (1 - \theta_i)^{n_i - y_i}\] である。

したがって、事後分布は、 \[p(\alpha, \beta | y, n, x) \propto p(\alpha, \beta | n, x)p(y | \alpha, \beta, n, x) \propto p(\alpha, \beta) \prod_{i=1}^k p(y_i | \alpha, \beta, n_i, x_i)\] となる。

\(\alpha\) と \(\beta\) についての事前情報がないので、無情報事前分布を使うことにする。 \(p(\alpha, \beta) \propto 1\) とすると、事後分布は、 \[p(\alpha, \beta | y, n, x) \propto \prod_{i=1}^k p(y_i | \alpha, \beta, n_i, x_i)\] となる。

大雑把な推定

最尤法を使って大雑把に推定を行うと、以下のような結果が得られる。

fit.logit <- glm(cbind(y, n - y) ~ x, data = Bio,
                 family = binomial(link = "logit"))
summary(fit.logit)
## 
## Call:
## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"), 
##     data = Bio)
## 
## Deviance Residuals: 
##        1         2         3         4  
## -0.17236   0.08133  -0.05869   0.12237  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.8466     1.0191   0.831    0.406
## x             7.7488     4.8728   1.590    0.112
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 15.791412  on 3  degrees of freedom
## Residual deviance:  0.054742  on 2  degrees of freedom
## AIC: 7.9648
## 
## Number of Fisher Scoring iterations: 7

等高線図を描く

\(\alpha\) と \(\beta\) の格子を作り、等高線図を描いてみよう。 ただし、\((\alpha, \beta) \in [-5, 10] \times [-10, 40]\) とする。

まず、格子を作る。

grd.wd.a = .1
grd.wd.b = .25
alpha.grid <- seq(-5, 10, by = grd.wd.a)
beta.grid <- seq(-10, 40, by = grd.wd.b)
Grid <- expand.grid(alpha.grid, beta.grid)
names(Grid) <- c("alpha", "beta")

次に、格子の各点上の事後確率を計算する。

invlogit <- function(x){ ## ロジスティック関数を定義する
    ## 引数：x = 数値ベクトル
    ## 返り値：y
    y <- 1 / (1 + exp(-x))
    return(y)
}

bio_post <- function(grid, data) { ## 式 (3.16) ただし、p(alpha, beta) = 1
    alpha <- grid$alpha
    beta <- grid$beta
    n <- length(alpha)
    ret_vec <- rep(NA, n)
    for (i in 1:n) {
        lik <- (invlogit(alpha[i] + beta[i] * data$x))^data$y *
              (1 - invlogit(alpha[i] + beta[i] * data$x))^(data$n - data$y)
        ret_vec[i] <- prod(lik)
    }
    return(ret_vec)
}
## unnormalized posterior density
Grid$density <- bio_post(grid = Grid, data = Bio)

得られた事後確率（ただし、まだ正規化されていない）を等高線図に描く。

bio.cont <- ggplot(Grid, aes(x = alpha, y = beta, z = density)) + geom_contour()
bio.cont + xlim(-5, 10) + ylim(-10, 40) +
    ggtitle("Contour plot of the posterior alpha and beta")

3次元図にしてみる。（web上には表示されていないので注意。）

plot3d(Grid$alpha, Grid$beta, Grid$density)

同時事後分布からのサンプリング

格子を利用し、確率を正規化する。

まず、\(\alpha\) の周辺事後分布を求める。

post.mg.alpha <- ddply(Grid, .(alpha), summarize, density = sum(density))
post.mg.alpha$density <- with(post.mg.alpha, density / sum(density))
bio.post.alpha <- ggplot(post.mg.alpha, aes(x = alpha, y = density)) + 
    geom_histogram(stat = "identity")
bio.post.alpha + ggtitle("Posterior marginal of alpha")

次に、事後分布から、\(\alpha\)、\(\beta\) を抽出する。 抽出は、まず \(p(\alpha |y)\) から\(\alpha\) を引き、抽出した \(\alpha\) をつかって \(p(\beta | \alpha, \beta)\) を引くという作業を行う。 その後、格子の隙間の値も取れるように、\(\alpha\)、\(\beta\) を少しだけずらす (jitter する)。

Post <- data.frame(alpha = rep(NA, 1000),
                   beta = rep(NA, 1000))
set.seed(2015-05-20)
for (s in 1:1000) {
    a <- sample(post.mg.alpha$alpha, size = 1,
                prob = post.mg.alpha$density)
    cond.Grid <- subset(Grid, alpha == a)
    p <- with(cond.Grid, density / sum(density))
    b <- sample(cond.Grid$beta, size = 1, prob = p)
    Post$alpha[s] <- a + runif(1, min = - grd.wd.a / 2, max = grd.wd.a / 2)
    Post$beta[s] <- b + runif(1, min = - grd.wd.b / 2, max = grd.wd.b / 2)
}

抽出した1000個の \((\alpha, \beta)\) を散布図にしてみよう。

bio.post <- ggplot(Post, aes(x = alpha, y = beta)) + geom_point()
bio.post + xlim(-5, 10) + ylim(-10, 40) +
    ggtitle("1000 draws from the posterior")

これを、上で描いた等高線図と比較せよ。

LD50 の事後分布

被験体の50%が死亡する確率 LD50 の事後分布を知りたいとする。 死亡確率が50%、つまり\(\theta_i = 0.5\) だから、 \[\mathrm{logit}^{-1}(\alpha + \beta x_i) = 0.5\] である。ロジットの逆関数が0.5 になるのは、\(\alpha + \beta x_i = 0\) のとき、すなわち、\(x_i = -\alpha / \beta\) である。したがって、私たちが知りたいのは、\(-\alpha / \beta\) の事後分布である。

ただし、薬物が死亡確率を上げない場合、死亡率を50%にする薬物量を考える意味がわからないので、\(\beta > 0\) の場合のみを考える。

この事後分布は、以下のようになる。

bio.ld50 <- ggplot(subset(Post, beta > 0), aes(x = -(alpha/beta))) +
    geom_histogram(aes(y = ..density..))
bio.ld50 + ggtitle("Draws from the posterior of the LD50") + xlab("LD50")



授業の内容に戻る