MCMC再訪(5): パス解析：Taglibro de H：SSブログ

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

MCMC再訪(5): パス解析　[統計] [編集]

前回の続き。今回はパス解析っぽいモデル。図示すると下の図のようになる。ただし、Zはポアソン分布。

解説は後で足すかも。

sample4.csv

"x","y","z"
"11.15","18.16","1"
"10.18","13.64","5"
"11.39","12.37","8"
" 3.91"," 7.49","7"
" 9.72"," 9.67","6"
" 7.17"," 8.63","11"
"10.16","13.84","4"
" 9.61","13.66","0"
"11.41","12.49","6"
"12.42","14.54","11"
"11.29","16.50","1"
" 7.20","13.85","5"
" 9.37","12.40","6"
"10.79","14.72","4"
" 9.99","14.79","2"
" 8.66","13.26","2"
"13.58","18.83","2"
" 8.72","13.03","0"
" 8.07","13.81","5"
" 9.11","13.99","3"
"10.03","14.12","2"
"13.95","15.38","3"
"10.29","11.94","7"
" 9.33","13.52","4"
"12.04","16.63","2"
"10.24","15.68","6"
"13.85","15.71","3"
"11.99","13.55","9"
" 8.67","13.31","2"
" 9.21","13.80","2"
" 8.88","13.57","2"
"10.07","13.53","8"
"11.13","14.55","6"
" 7.42"," 9.98","6"
"11.34","16.53","3"
" 7.68","10.07","3"
"12.25","17.32","2"
"10.94","14.19","3"
"10.09","14.40","3"
"11.04","15.14","2"
" 8.07","11.23","6"
"10.00","15.91","1"
"10.67","14.54","8"
"12.32","17.44","0"
"12.59","15.37","2"
"11.85","17.12","0"
"12.61","16.66","4"
"11.66","14.33","1"
"14.74","16.86","8"
"10.58","16.18","2"

sample4.bug

var
  N,
  X[N],
  Y[N],
  Z[N],
  Y.hat[N], mu.z[N],
  X.bar[N], Y.bar[N],
  beta.y, beta.yx,
  beta.z, beta.zx, beta.zy,
  tau.y, sigma.y;

model {
  for (i in 1:N) {
    Z[i] ~ dpois(mu.z[i]);
    log(mu.z[i]) <- beta.z + beta.zx * (X[i] - X.bar) +
                             beta.zy * (Y[i] - Y.bar);
    Y[i] ~ dnorm(Y.hat[i], tau.y);
    Y.hat[i] <- beta.y + beta.yx * (X[i] - X.bar);
  }

  # centering
  X.bar <- mean(X[]);
  Y.bar <- mean(Y[]);

  # priors
  beta.y ~ dnorm(0, 1.0E-6);
  beta.z ~ dnorm(0, 1.0E-6);
  beta.yx ~ dnorm(0, 1.0E-6);
  beta.zx ~ dnorm(0, 1.0E-6);
  beta.zy ~ dnorm(0, 1.0E-6);
  tau.y ~ dgamma(1.0E-3, 1.0E-3);

  sigma.y <- 1/sqrt(tau.y);
  ief.zx <- beta.zy * beta.yx;
  tef.zx <- beta.zx + beta.zy * beta.yx;
}

sample4.R

##
## Sample 4
##

## WinBUGS
library(R2WinBUGS)

## 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
}

## data
data <- read.csv("sample4.csv")

## number of chains
n.chains <- 3

## initial values
init1 <- list(beta.y =  1, beta.z = 10, beta.yx = 2,
              beta.zx = 1, beta.zy = 1, tau.y = 1)
init2 <- list(beta.y =  5, beta.z =  1, beta.yx = 4,
              beta.zx = 5, beta.zy = 5, tau.y = 0.1)
init3 <- list(beta.y =  10, beta.z = 10, beta.yx = 0,
              beta.zx = 0, beta.zy = 0, tau.y = 10)

## MCMC
post.samp <- bugs(data = list(N = N, X = x, Y = y, Z = z),
                  inits = list(init1, init2, init3),
                  parameters.to.save = c("beta.y", "beta.z",
                                         "beta.yx",
                                         "beta.zx", "beta.zy",
                                         "sigma.y",
                                         "ief.zx", "tef.zx"),
                  model = "sample4.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)

## show results
print(post.samp, digits = 3)