※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 04月 11日
N <- 20 linkFunction <- function(b1, b2) { function(x) { 1 / (1 + exp(-b1 - b2 * x)) } } distribution <- function(x, linkFunction) { l <- match.fun(linkFunction) function(y) { dbinom(y, N, l(x)) } } xs <- seq(0, 10, by = 2) ys <- 0:N xl <- "y" yl <- "Probability" legends <- paste0("x = ", xs) for (b1 in seq(-0.5, 0.5, by = 0.2)) { for (b2 in seq(-0.5, 0.5, by = 0.2)) { l <- linkFunction(b1, b2) title <- paste0("q = 1 / (1 + exp(-(", b1, " + ", b2, "x)))") plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl) for (x in xs) { d <- distribution(x, l) lines(ys, d(ys), type = "l") points(ys, d(ys), pch = x) } legend("topright", legend = legends, pch = xs) } }リンク関数はロジスティック関数です。 コードでは以下の部分です。
linkFunction <- function(b1, b2) { function(x) { 1 / (1 + exp(-b1 - b2 * x)) } }これが二項分布のパラメータの一つである事象の発生確率 `q` を与えます。 もう一つのパラメータは試行回数 `N` ですが、これは `20` にしました。 ということで、リンク関数のパラメータ `b1` と `b2` さえ与えれば、説明変数 `x` に対して応答変数 `y` の分布がどのように変化するかが分かります。 コードを実行すると様々な `b1` と `b2` の値に対してモデルをプロットします。 例えば下図のような図がプロットされます。