r - R2OpenBUGS: 中間結果からの最終治療結果のモデル化

Question

R2OpenBUGS を使っている人はいますか? r2winbugs を使用する必要がありますか? ...

私は、(単一の) 中間 (3 か月) の転帰を持つ患者のサンプルに対して、最終 (2 年間) の治療転帰 (成功、死亡、デフォルトまたは失敗など) をモデル化しようとしています。

R2OpenBUGS は多項式ノードで奇妙な事後分布を示しています。結果のうち 2 つは一定で、他の 2 つの結果は等しく、結果の総数はコホート サイズよりも大きくなっています。

私は何を間違っていますか？よろしくお願いします！コードアウトの出力は以下のとおりです。

   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)

出力の一部を次に示します。

> 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).