MCMC再訪(2): ロジスティック回帰：Taglibro de H：SSブログ

	ブログをはじめるログイン

MCMC再訪(2): ロジスティック回帰　[統計] [編集]

(1)につづき、今回はロジスティック回帰をやってみる。

ある生物の集団があったとする。K = 10個体ずつに、X (X₁ = 0, X₂ = 1, X₃ = 2,..., X₁₁ = 10)単位だけある薬品を与えたところ、それぞれについて10個体中のY_i個体が死亡したとする。このとき、XとYとの関係をモデル化し、パラメータを推定する。

以下のようなコード

library(R2WinBUGS)

R2WinBUGSパッケージを使用する。

## data
K <- 10
X <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
Y <- c(1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10)
N <- length(X)

データを用意。

model <- function() {
  for (i in 1:N) {
    Y[i] ~ dbin(p[i], K);
    logit(p[i]) <- beta + beta.x * X[i];
  }
  
  # priors
  beta ~ dnorm(1.0E-6, 1.0E-6);
  beta.x ~ dnorm(1.0E-6, 1.0E-6);
}
write.model(model, "sample2.bug")

モデル。Y[i]は、試行数K、確率p[i]の二項分布に従うとする。p[i]とX[i]との関係が次の行。betaおよびbeta.xはパラメータで、事前分布は無情報とする。

## Settings
platform <- .Platform$OS.type
if (platform == "unix") {
  useWINE = TRUE
  ## check Darwine
  darwine.bin <- "/Applications/Darwine/Wine.bundle/Contents/bin"
  if (file.exists(darwine.bin)) {
    wine <- paste(darwine.bin, "wine", sep = "/")
    winepath <- paste(darwine.bin, "winepath", sep = "/")
  } else {
    ## use `which'
    wine <- system("which wine", intern = TRUE)
    winepath <- system("which winepath", intern = TRUE)
    if (wine == ""|| winepath == "") {
      stop("Wine not found.")
    }
  }
  bugs.dir <- paste(Sys.getenv("HOME"),
                    ".wine/drive_c/Program Files/WinBUGS14",
                    sep = "/")
  if (!file.exists(bugs.dir)) {
    stop("bugs directory not found.")
  }
} else if (platform == "windows") {
  ## Windows
  useWINE = FALSE
  bugs.dir <- "C:/Program Files/WinBUGS14"
  if (!file.exists(bugs.dir)) {
    stop("bugs directory not found.")
  }
  wine <- NULL
  winepath <- NULL
}

実行環境に応じた設定をおこなう。

## number of chains
n.chains <- 3

## initial values
init1 <- list(beta = -10, beta.x = 0)
init2 <- list(beta =  -5, beta.x = 2)
init3 <- list(beta =   0, beta.x = 4)

chainの数と、初期値を用意する。

## MCMC
post.samp <- bugs(data = list(X = X, Y = Y, N = N, K = K),
                  inits = list(init1, init2, init3),
                  parameters.to.save = c("beta", "beta.x"),
                  model = "sample2.bug",
                  n.chains = n.chains,
                  n.iter = 21000,
                  n.burnin = 1000,
                  n.thin = 10,
                  debug = FALSE,
                  bugs.directory = bugs.dir,
                  working.directory = getwd(),
                  clearWD = TRUE,
                  useWINE = useWINE,
                  newWINE = TRUE,
                  WINE = wine,
                  WINEPATH = winepath)

WinBUGSの実行

print(post.samp, digits = 3)

結果の表示。以下の結果が表示される。

Inference for Bugs model at "sample2.bug", fit using WinBUGS,
 3 chains, each with 21000 iterations (first 1000 discarded), n.thin = 10
 n.sims = 6000 iterations saved
           mean    sd   2.5%    25%    50%    75%  97.5%  Rhat n.eff
beta     -2.161 0.492 -3.164 -2.484 -2.144 -1.823 -1.246 1.001  3000
beta.x    0.556 0.101  0.372  0.486  0.552  0.621  0.763 1.001  6000
deviance 30.078 2.024 28.140 28.670 29.460 30.880 35.481 1.001  6000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.0 and DIC = 32.0
DIC is an estimate of expected predictive error (lower deviance is better).

## mcmc object
post.mcmc <- mcmc.list(lapply(1:(post.samp$n.chain),
                              function(chain)
                              as.mcmc(post.samp$sims.array[,chain,])))
plot(post.mcmc)

軌跡と密度を表示する。

beta <- post.samp$mean$beta
beta.x <- post.samp$mean$beta.x

x <- seq(0, 10, len = 100)
logit.p <- beta + beta.x * x
exp.Y <- K * exp(logit.p) / (1 + exp(logit.p))
plot(X, Y, type = "p", ylim = c(0, 10), las = 1)
lines(x, exp.Y)

データと、MCMCから求められた期待とを重ねてプロットする。