[R] Getting WinBUGS Leuk example to work from R using R2winBUGS

[Uwe Ligges]

On 17.02.2013 08:47, Andy Cox wrote:
> I am trying to learn to use winBUGS from R, I have experience with R.
> I have managed to successfully run a simple example from R with no
> problems. I have been trying to run the Leuk: Survival from winBUGS
> examples Volume 1. I have managed to run this from winBUGS GUI with no
> problems. My problem is try as I might ( and I have been trying and
> searching for days) I cannot get it to run using R2winBUGS.I am sure
> it is something simple.
>
> The error message I get if I try and set inits in the script is
>
> Error in bugs(data = L, inits = inits,
>        parameters.to.save = params, model.file "model.txt",  :
>        Number of initialized chains (length(inits)) != n.chains
>
> I know this means I have not initialised some of the chains, but I am
> pasting the inits code from winbugs examples manual and all other
> settings seem to me to be the same as when run on the winBUGS GUI.
>
>
>
> If I try inits=NULL I get another error message
>
>      display(log)
> check(C:/BUGS/model.txt)
> model is syntactically correct
> data(C:/BUGS/data.txt)
> data loaded
> compile(1)
> model compiled
> gen.inits()
> shape parameter (r) of gamma dL0[1] too small -- cannot sample
> thin.updater(1)
> update(500)
> command #Bugs:update cannot be executed (is greyed out)
> set(beta)
>
> Which indicates to me I will still have problems after solving the
> first one!! I am about to give up on using winBUGS, please can someone
> save me? I know I am probably going to look stupid, but everyone has
> to learn:-)
>
> I have also tried changing nc<-2 (On advice, which doesnt work and
> gives an uninitialised chain error)
>
> I am using winBUGS 1.4.3 on Windows XP 2002 SP3
>
> My R code is below, many thanks for at least reading this far.
>
>     rm(list = ls())
>
> L<-list(N = 42, T = 17, eps = 1.0E-10,
>          obs.t = c(1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12,
> 12, 15, 17, 22, 23, 6,
>          6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32,
> 32, 34, 35),
>          fail = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0,
> 0),
>          Z = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
>                 -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
> -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
> -0.5),
>          t = c(1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35))
>
> ### 5.4. Analysis using WinBUGS
> library(R2WinBUGS)      # Load the R2WinBUGS library CHOOSE to use WinBUGS
>      #library(R2OpenBUGS)        # Load the R2OpenBUGS library CHOOSE
> to use OpenBUGS
> setwd("C://BUGS")
>
> # Save BUGS description of the model to working directory
> sink("model.txt")
> cat("
>
> model
> {
> # Set up data
> for(i in 1:N) {
> for(j in 1:T) {
>
> # risk set = 1 if obs.t >= t
> Y[i,j] <- step(obs.t[i] - t[j] + eps)
> # counting process jump = 1 if obs.t in [ t[j], t[j+1] )
> # i.e. if t[j] <= obs.t < t[j+1]
> dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i]
> }
> }
> # Model
> for(j in 1:T) {
> for(i in 1:N) {
> dN[i, j] ~ dpois(Idt[i, j]) # Likelihood
> Idt[i, j] <- Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity
> }
> dL0[j] ~ dgamma(mu[j], c)
> mu[j] <- dL0.star[j] * c # prior mean hazard
> # Survivor function = exp(-Integral{l0(u)du})^exp(beta*z)
> S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5));
> S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5));
> }
> c <- 0.001
> r <- 0.1
> for (j in 1 : T) {
> dL0.star[j] <- r * (t[j + 1] - t[j])
> }
> beta ~ dnorm(0.0,0.000001)
> }
>
>
> ",fill=TRUE)
> sink()
>
> params<- c("beta","S.placebo","S.treat")
>
> inits<-list( beta = 0.0,
>       dL0 = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
>               1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0))


You need a list containing one list of initial values for each chain. 
For you nc=1 number of chains, this means:

  inits <- list(list( beta = 0.0,
        dL0 = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0)))


Best,
Uwe Ligges


> # MCMC settings
> nc <-1                  # Number of chains
> ni <- 1000              # Number of draws from posterior (for each chain)
> ns<-1000 #Number of sims (n.sims)
> nb <- floor(ni/2)                   # Number of draws to discard as burn-in
> nt <- max(1, floor(nc * (ni - nb) / ns))# Thinning rate
> Lout<-list()
>
>
>
> # Start Gibbs sampler: Run model in WinBUGS and save results in object
> called out
> out <- bugs(data = L, inits =inits, parameters.to.save = params,
> model.file = "model.txt",
> n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = T, DIC
> = TRUE,digits=5,
> codaPkg=FALSE, working.directory = getwd())
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>