R2OpenBUGS: modelling final treatment outcomes from intermediate outcomes

Question

Is anyone using R2OpenBUGS? Should I rather be using r2winbugs? ...

I am trying to model final (2-year) treatment outcomes (e.g. success, death, default or failure) for my sample of patients with a (single) intermediate (3-month) outcome.

R2OpenBUGS is giving me a strange posterior on the multinomial node, in which two of the outcomes are constant, the other two outcomes are equal, and the total number of outcomes is greater than the cohort size.

What am I doing wrong? Many thanks in advance! Code out output are below.

   library(R2OpenBUGS)

    model <- function() { 
      # Prior : distribution of final outcomes for treatment cohort N_tx
      outc[1:4] ~ dmulti(p.outc[],N_tx)
      p.outc[1] <- 164/1369
      p.outc[2] <- 907/1369
      p.outc[3] <- 190/1369
      p.outc[4] <- 108/1369
      # Prior : distribution of intermediate outcome (proxy of final outcome) for each final outcome cohort 
      # (e.g. proportion of patient with final outcome 1 that exhibited the intermediate outcome)
      cr_1 ~ dunif(0.451, 0.609)
      cr_2 ~ dunif(0.730, 0.787)
      cr_3 ~ dunif(0.559, 0.700)
      cr_4 ~ dunif(0.148, 0.312)
      # Probability p of the intermediate outcome given prior distributions
      p <- (outc[1]*cr_1+outc[2]*cr_2+outc[3]*cr_3+outc[4]*cr_4)/N_tx
      # Likelihood function for the number of culture conversions at 3 months among those still on treatment in month 6 (excludes confirmed deaths and defaulters)
      cs ~ dbin(p,N_tx)
    } 

    # N_tx is the number of patients in our cohort 
    N_tx <- 100
    # cs is the number of patient exhibiting the intermediate outcome
    cs <- 80

    data <- list("N_tx", "cs")

    inits <- function() { list(outc=c(round(164/1369*N_tx),
                                      round(907/1369*N_tx),
                                      round(190/1369*N_tx),
                                      round(108/1369*N_tx)),
                               cr_1=87/(87+77), 
                               cr_2=689/(689+218), 
                               cr_3=120/(120+70), 
                               cr_4=24/(24+84))
    }

    params <- c("outc") 

    model.file <- file.path(tempdir(), "model.txt") 
    write.model(model, model.file)

    out <- bugs(data, inits, params, model.file, n.iter=100000,debug=TRUE)

    all(out$summary[,"Rhat"] < 1.1) 

    out$mean["outc"] 
    out$sd["outc"] 

    print(out, digits=5)

And here are some of the outputs:

> all(out$summary[,"Rhat"] < 1.1) 
[1] TRUE
> 
> out$mean["outc"] 
$outc
[1] 15.53095 66.00000 14.00000 15.53095

> out$sd["outc"] 
$outc
[1] 3.137715 0.000000 0.000000 3.137715

> 
> print(out, digits=5)
Inference for Bugs model at "C:\", 
Current: 3 chains, each with 1e+05 iterations (first 50000 discarded)
Cumulative: n.sims = 150000 iterations saved
             mean      sd     2.5%    25%    50%    75% 97.5%    Rhat  n.eff
outc[1]  15.53095 3.13771 10.00000 13.000 15.000 18.000 22.00 1.00100 130000
outc[2]  66.00000 0.00000 66.00000 66.000 66.000 66.000 66.00 1.00000      1
outc[3]  14.00000 0.00000 14.00000 14.000 14.000 14.000 14.00 1.00000      1
outc[4]  15.53095 3.13771 10.00000 13.000 15.000 18.000 22.00 1.00100 130000
deviance  8.59096 2.23382  5.08097  6.927  8.323  9.963 13.66 1.00102  55000

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 = var(deviance)/2)
pD = 2.5 and DIC = 11.1
DIC is an estimate of expected predictive error (lower deviance is better).