※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 05月 30日
N <- 20 linkFunction <- function(b1, b2, x) { function(r) { 1 / (1 + exp(-b1 - b2 * x - r)) } } distribution <- function(y, linkFunction) { l <- match.fun(linkFunction) function(r) { dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = 2) } } xl <- "r" yl <- "Probability" ys <- seq(0, 20, by = 5) rs <- seq(-5, 5, by = 0.5) legends = paste0("y = ", ys) for (b1 in seq(-0.5, 0.5, by = 0.2)) { for (b2 in seq(-0.5, 0.5, by = 0.2)) { for (x in 0:3) { l <- linkFunction(b1, b2, x) title <- paste0("logit(q) = ", b1, " + ", b2, " * ", x, " + r, r ~ N(0, 2)") plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.1), main = title, xlab = xl, ylab = yl) for (i in 1:5) { d <- distribution(ys[i], l) lines(rs, d(rs), type = "l") points(rs, d(rs), pch = ys[i]) } legend("topright", legend = legends, pch = ys) } } }リンク関数は `logit(q) = b1 + b2 * x + r` です。 `r` は平均 `0`、標準偏差 `2` の正規分布から生成される変数です。 事象の試行回数は `20` にしました。 リンク関数を使って二項分布する事象の発生確率 `q` が求まれば、事象が `y` 回起こる確率を計算できます。 あとは、正規分布から `r` が生成される確率と、二項分布を掛け合わせれば、欲しい値になります。 様々な `b1`、`b2`、`x` を与えて、`r` に対して確率をプロットしてみました。 コードを実行すると144枚のグラフがプロットされます。 例えば次のようなグラフが得られます。