※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 05月 16日
N <- 20 linkFunction <- function(b1, b2) { function(x, r) { 1 / (1 + exp(-b1 - b2 * x - r)) } } distribution <- function(x, r, linkFunction) { l <- match.fun(linkFunction) function(y) { dbinom(y, N, l(x, r)) } } rs <- -2:2 ys <- 0:N pchs <- 1:5 xl <- "y" yl <- "Probability" legends <- paste0("r = ", rs) for (b1 in c(-0.5, 0.5, by = 0.2)) { for (b2 in c(-0.5, 0.5, by = 0.2)) { l <- linkFunction(b1, b2) title <- paste0("logit(q) = ", b1, " + ",b2," * 2"," + r") plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl) for (i in 1:5) { d <- distribution(2, rs[i], l) lines(ys, d(ys), type = "l") points(ys, d(ys), pch = pchs[i]) } legend("topright", legend = legends, pch = pchs) } }このモデルでは、二項分布のパラメータである事象の発生確率 `q` が以下のコードで与えられます。
linkFunction <- function(b1, b2) { function(x, r) { 1 / (1 + exp(-b1 - b2 * x - r)) } }GLMのモデルと比べると、新しい変数として `r` が追加されています。 この `r` が説明変数 `x` とは違った、何とは特定されていない、余分な事象の発生確率への影響を表します。 コードを実行すると様々な `b1` と `b2` と `r` の値に対してモデルをプロットします。 `x` は `2` で固定にしました。 事象の試行回数は `N = 20` で、事象の発生回数は `y` です。 例えば下図のような図がプロットされます。