R: runjags：Taglibro de H：SSブログ

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

R: runjags　[統計] [編集]

runjags: Run Bayesian MCMC Models in the BUGS syntax from Within R

JAGSを呼び出すパッケージ。収束状況をみて自動的に繰り返し回数を決定する関数もある。

ためしてみた。

library(MASS)
data(bacteria)

n <- nrow(bacteria)
yInt <- as.integer(bacteria$y == "y")
week <- bacteria$week
trtdrug <- as.integer(bacteria$trt == "drug")
trtdrugplus <- as.integer(bacteria$trt == "drug+")
id <- as.integer(bacteria$ID)
n.id <- length(levels(bacteria$ID))

model <- "
model {
  for (i in 1:n) {
    yInt[i] ~ dbern(p[i])
    logit(p[i]) <- b + b.week * week[i] +
                   b.drug * trtdrug[i] +
                   b.drugplus * trtdrugplus[i] +
                   e[id[i]]
  }
  for (j in 1:n.id) {
    e[j] ~ dnorm(0.0, tau)
  }
  b ~ dnorm(0.0, 1.0E-4)
  b.week ~ dnorm(0.0, 1.0E-4)
  b.drug ~ dnorm(0.0, 1.0E-4)
  b.drugplus ~ dnorm(0.0, 1.0E-4)
  tau ~ dgamma(1.0E-3, 1.0E-3)
  sigma <- 1/sqrt(tau)
}
"

vars <- c("b", "b.week", "b.drug", "b.drugplus", "sigma")
data <- list(yInt = yInt, week = week,
             trtdrug = trtdrug, trtdrugplus = trtdrugplus,
             id = id, n = n, n.id = n.id)
init1 <- list(.RNG.seed = 123,
              .RNG.name = "base::Mersenne-Twister",
              b = -10, b.week = -3, b.drug = 0,
              b.drugplus = 10, tau = 1)
init2 <- list(.RNG.seed = 1234,
              .RNG.name = "base::Mersenne-Twister",
              b = -3, b.week = 0, b.drug = 3,
              b.drugplus = 5, tau = 10)
init3 <- list(.RNG.seed = 12345,
              .RNG.name = "base::Mersenne-Twister",
              b = 10, b.week = -5, b.drug = 10,
              b.drugplus = 10, tau = 0.1)
inits <- list(init1, init2, init3)


library(runjags)
results <- autorun.jags(model = model, monitor = vars,
                        data = data, inits = inits,
                        n.chains = 3)

結果

Auto-run JAGS

Running a pilot chain...
Calling the simulation... (this may take some time)
Welcome to JAGS 2.1.0 on Fri Jul  2 21:11:47 2010
JAGS is free software and comes with ABSOLUTELY NO WARRANTY
Loading module: basemod
Loading module: bugs
. <>. <Reading data file data.txt
>. <>Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 2923
. <Reading parameter file inits1.txt
><>. <Reading parameter file inits2.txt
><>. <Reading parameter file inits3.txt
><>. . <>Updating 5000
---------------------------------------| 5000
**************************************** 100%
. <>. <>. <>. <>. <>. <>Updating 10000
---------------------------------------| 10000
**************************************** 100%
. <>. <><>. <><>. <><>. 
Simulation complete.  Reading coda files...
Coda files loaded successfully

Calculating the Gelman-Rubin statistic....
The Gelman-Rubin statistic is below 1.05 for all parameters

Calculating the necessary sample length based on the Raftery and Lewis's diagnostic...
IMPORTANT:  The sample size(s) of monitored node(s) 'b' & 'b.week' & 'b.drug' & 'b.drugplus' & 'sigma' have a high autocorrelation dependance in chain(s) 1 & 2 & 3.  Re-running the model with a different formulation or better initial values may help to reduce autocorrelation.
The model will need to be run for a further 64404 updates.  This will take approximately 3.1 minutes.
Continue with the simulation?  y

Calling the simulation... (this may take some time)
Welcome to JAGS 2.1.0 on Fri Jul  2 21:12:55 2010
JAGS is free software and comes with ABSOLUTELY NO WARRANTY
Loading module: basemod
Loading module: bugs
. <>. <Reading data file data.txt
>. <>Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 2923
. <Reading parameter file inits1.txt
><>. <Reading parameter file inits2.txt
><>. <Reading parameter file inits3.txt
><>. . <>Updating 200
---------------------------------------| 200
**************************************** 100%
. <>. <>. <>. <>. <>. <>Updating 64404
---------------------------------------| 64400
**************************************** 100%
* 100%
. <>. <><>. <><>. <><>. 
Simulation complete.  Reading coda files...
Coda files loaded successfully
Necessary sample length achieved
REMINDER:  There was a high degree of autocorrelation for 2 parameter(s)
Auto-run JAGS complete.

*PLEASE NOTE:  THIS SOFTWARE IS INTENDED FOR EDUCATIONAL PURPOSES ONLY*
*YOU SHOULD ASSESS CONVERGENCE AND AUTOCORRELATION MANUALLY BEFORE RELYING ON RESULTS PROVIDED*

>  results$summary

Iterations = 1:74404
Thinning interval = 7 
Number of chains = 3 
Sample size per chain = 10630 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean     SD  Naive SE Time-series SE
b           3.4257 0.7365 0.0041240      0.0117294
b.week     -0.1552 0.0543 0.0003041      0.0005316
b.drug     -1.4346 0.7727 0.0043269      0.0085372
b.drugplus -0.8984 0.7766 0.0043490      0.0081712
sigma       1.4218 0.5288 0.0029610      0.0122668

2. Quantiles for each variable:

              2.5%     25%     50%     75%    97.5%
b           2.1706  2.9036  3.3597  3.8800  5.05778
b.week     -0.2650 -0.1906 -0.1540 -0.1176 -0.05277
b.drug     -3.0815 -1.8988 -1.4018 -0.9293  0.01208
b.drugplus -2.5305 -1.3751 -0.8691 -0.3945  0.56706
sigma       0.3557  1.0845  1.3966  1.7426  2.53068

>  results$autocorr
                  b       b.week      b.drug  b.drugplus      sigma
Lag 0    1.00000000 1.0000000000  1.00000000  1.00000000 1.00000000
Lag 7    0.69562980 0.2426710270  0.46920553  0.43361715 0.72048587
Lag 35   0.26084934 0.0612951207  0.09253383  0.09093770 0.34937331
Lag 70   0.09835150 0.0303896344  0.03303743  0.02875206 0.19673322
Lag 350 -0.00639729 0.0006886601 -0.00803727 -0.01733838 0.02834677
>  results$psrf
Potential scale reduction factors:

           Point est. 97.5% quantile
b                1.00           1.00
b.week           1.00           1.00
b.drug           1.00           1.00
b.drugplus       1.00           1.00
sigma            1.00           1.01

Multivariate psrf

1.00

Target psrf

1.05