0

在尝试使用 调整工作WinBUGS模型并将其缓解为 RR2WinBUGS时，我收到了几条错误消息。我相信与 的规范有关inits，但一直无法解决问题。

第一个是错误消息正在使用inits1：

list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0))

Error in bugs(mydata, inits = inits1, model.file = "mtcfe.txt", parameters = c("or"),  : 
Number of initialized chains (length(inits)) != n.chains

在阅读了 Uwe Liggers “列表包含列表解决方案”的建议修复后，我将 inits1 修改为inits2：

inits2 <- list(list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0)))

并收到错误：

Error in bugs.run(n.burnin, bugs.directory, WINE = WINE, useWINE = useWINE,  :  Look at the log file and try again with 'debug=TRUE' to figure out what went wrong within Bugs.

我还尝试了 AndyC 在这篇文章“getting-winbugs-leuk-example-to-work-from-r-using-r2winbugs”中建议的修复。通过更改inits为：

inits4 <- function(){list(list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0)))}

并收到错误：

Error in bugs.run(n.burnin, bugs.directory, WINE = WINE, useWINE = useWINE,  :  Look at the log file and try again with 'debug=TRUE' to figure out what went wrong within Bugs.

这是我对R2WinBUGS. 请注意，我inits在代码中包含了几个不起作用的代码：

work.dir <- "removed from example"
setwd(work.dir)
getwd()      # check working directory

# Load Package
library(R2WinBUGS)

# Read data
mydata <- list(nt=5,
     ns=4,
          r=structure(
       .Data = c(2506,7834,6729,2139,
                 2548,7860,6710,4418),
       .Dim = c(4,2)),
     n=structure(
       .Data = c(2697, 8212, 7266, 2333,
                 2701,8280,7257,4687),
       .Dim= c(4,2)),
     t=structure(
       .Data = c( 1,1,1,1,
                  2,3,4,5),
       .Dim = c(4,2)),
     na=structure(
       .Data = c(2,2,2,2))
        )

bugs.data(mydata)

# Set initial values

inits <- function(){list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0))}
#inits2 <- list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0))
#inits3 <- function(){list(inits2,inits2)}
#inits4 <- function(){list(inits,inits)} 
#inits5 <- list(inits,inits)
# CALL WinBUGS AND SAVE RESULTS IN VARIABLE out.re
out.re <- bugs(mydata,inits=inits1,         # load data and initial values
               model.file="mtcfe.txt",      # file with model to run
               parameters=c("or"),
               n.thin=1,
               n.chains=1, n.iter=1500, n.burnin=500,
               bugs.directory=bd, 
               working.directory=work.dir,
               debug=TRUE)

print(out.re,digits=4)   #lists all results as in WinBUGS stats(*)

这是工作WinBUGS代码，改编自WinBUGS布里斯托尔和莱斯特大学举办的课程：

# Binomial likelihood, logit link, MTC
# Fixed effect model
model{                                                                      # *** PROGRAM STARTS
  for(i in 1:ns){                                                               # LOOP THROUGH STUDIES
    mu[i] ~ dnorm(0,.0001)                                              # vague priors for all trial baselines
    for (k in 1:na[i]) {                                                       # LOOP THROUGH ARMS
      r[i,k] ~ dbin(p[i,k],n[i,k])                                           # binomial likelihood
      logit(p[i,k]) <- mu[i] + d[t[i,k]]-d[t[i,1]]                        # model for linear predictor
      rhat[i,k] <- p[i,k] * n[i,k]                                            # expected value of the numerators
      dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k]))              #Deviance contribution
                       + (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k])))
    }
    resdev[i] <- sum(dev[i,1:na[i]])                  # summed residual deviance contribution for this trial
  }
  totresdev <- sum(resdev[])                             #Total Residual Deviance
  d[1]<- 0                                                       # treatment effect is zero for reference treatment
  for (k in 2:nt)  { d[k] ~ dnorm(0,.0001) }           # vague priors for treatment effects

  # pairwise ORs and LORs for all possible pair-wise comparisons
  for (c in 1:(nt-1)) {  for (k in (c+1):nt) {
    or[c,k] <- exp(d[k] - d[c])
    lor[c,k] <- (d[k]-d[c])
  }
  }

  # ranking 
  for (k in 1:nt) {
    rk[k] <- nt+1-rank(d[],k)                      # assumes events are “good”
    # rk[k] <- rank(d[],k)                           # assumes events are “bad”
    best[k] <- equals(rk[k],1)                    #calculate probability that treat k is best
  }

  # Absolute effects
  A ~ dnorm(-2.6,precA)
  precA <- pow(0.38,-2)                           #  prior precision for Treatment A, sd=0.38 on logit scale
  for (k in 1:nt) { logit(T[k]) <- A + d[k] }
}                                                                                 # *** PROGRAM ENDS


#Inits
list(A=1, d=c(NA,0,0,0,0), mu=c(0,0,0,0))

#Data
list(nt=5.00000E+00, ns=4.00000E+00, r= structure(.Data= c(2.50600E+03, 2.54800E+03, 7.83400E+03, 7.86000E+03, 6.72900E+03, 6.71000E+03, 2.13900E+03, 4.41800E+03), .Dim=c(4, 2)), n= structure(.Data= c(2.69700E+03, 2.70100E+03, 8.21200E+03, 8.28000E+03, 7.26600E+03, 7.25700E+03, 2.33300E+03, 4.68700E+03), .Dim=c(4, 2)), t= structure(.Data= c(1.00000E+00, 2.00000E+00, 1.00000E+00, 3.00000E+00, 1.00000E+00, 4.00000E+00, 1.00000E+00, 5.00000E+00), .Dim=c(4, 2)), na=c(2.00000E+00, 2.00000E+00, 2.00000E+00, 2.00000E+00))
4

0 回答 0