在 R2WinBUGS 中调用 bugs() 函数失败日志中的错误消息说“使用了未定义的节点 n”

问题描述

模型和初始化直接从 NICE（TSD 2 2016 版，Dias 等人）的广泛使用的网络元分析银屑病示例 6（RE 条件二项式概率示例代码）复制而来，我正在尝试将其转换为使用在 R 中使用 R2WinBUGS 包。

我将示例中的研究数据传输到一个 .csv 文件中，并已读入以 R 格式发送到 WinBUGS。

代码，包括模型、数据和初始化设置。


library(R2WinBUGS)
#setwd
setwd("C:/Users/User1/Autotask Workplace/Axis Working Folder/Niamh Russell Folder/WinBUGS")
# read in csv file
rin<-read.csv("Example 6 TSU data.csv")

# Binomial likelihood,logit link,pairwise Meta-analysis (2 treatments)
# Random effects model
psoriasisRE.model<-function(){ # *** PROGRAM STARTS
  # Binomial likelihood,probit link (different categories)
  # Random effects model for multi-arm trials
  for(i in 1:ns){                      #   LOOP THROUGH STUDIES
    w[i,1] <- 0    # adjustment for multi-arm trials is zero for control arm
    delta[i,1] <- 0             # treatment effect is zero for control arm
    # vague priors for all trial baselines
    mu[i] ~ dnorm(0,.0001)
    for (k in 1:na[i]) {               # LOOP THROUGH ARMS
      p[i,k,1] <- 1                    # Pr(PASI >0)
      for (j in 1:nc[i]-1) {           # LOOP THROUGH CATEGORIES
        # binomial likelihood
        r[i,j] ~ dbin(q[i,j],n[i,j]) 
        # conditional probabilities
        q[i,j] <- 1-(p[i,C[i,j+1]]/p[i,j]])
        # linear predictor
        theta[i,j] <- mu[i] + delta[i,k] + z[C[i,j+1]-1]
        rhat[i,j] <-  q[i,j] * n[i,j]  # predicted number events 
        #Deviance contribution of each category
        dv[i,j] <- 2 * (r[i,j]*(log(r[i,j])-log(rhat[i,j]))  
                          +(n[i,j]-r[i,j])*(log(n[i,j]) - log(n[i,j]-rhat[i,j])))
      }
      dev[i,k] <- sum(dv[i,1:nc[i]-1]) # deviance contribution of each arm
      for (j in 2:nc[i])  {            # LOOP THROUGH CATEGORIES
        p[i,j]] <- 1 - phi.adj[i,j]   # link function
        # adjust link function phi(x) for extreme values that can give numerical errors
        # when x< -5,phi(x)=0,when x> 5,phi(x)=1: use only if needed
        phi.adj[i,j] <- step(5+theta[i,j-1])*(step(theta[i,j-1]-5)+ step(5-theta[i,j-1])*phi(theta[i,j-1]) )
      }
    }
    for (k in 2:na[i]) {               # LOOP THROUGH ARMS
      delta[i,k] ~ dnorm(md[i,k],taud[i,k])
      # mean of LHR distributions,with multi-arm trial correction
      md[i,k] <-  d[t[i,k]] - d[t[i,1]] + sw[i,k]
      # precision of LHR distributions (with multi-arm trial correction)
      taud[i,k] <- tau *2*(k-1)/k
      # adjustment,multi-arm RCTs
      w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]])
      # cumulative adjustment for multi-arm trials
      sw[i,k] <- sum(w[i,1:k-1])/(k-1)
    }
    #  summed residual deviance contribution for this trial
    resdev[i] <- sum(dev[i,1:na[i]])       
  }   
  z[1] <- 0                            # set z50=0
  # Set priors for z,for any number of categories
  for (j in 2:Cmax-1) {
    z.aux[j] ~ dunif(0,5)              # priors 
    z[j] <- z[j-1] + z.aux[j]          # ensures z[j]~Uniform(z[j-1],z[j-1]+5)
  }
  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
  sd ~ dunif(0,5)                      # vague prior for between-trial SD
  tau <- pow(sd,-2)     # between-trial precision = (1/between-trial variance)
  # Provide estimates of treatment effects T[k] on the natural 
  # (probability) scale
  # Given a Mean Effect,meanA,for 'standard' treatment A,# with precision (1/variance) precA
  A ~ dnorm(meanA,precA)
  for (k in 1:nt) { 
    # calculate prob of achieving PASI50,75,90 on treat k
    for (j in 1:Cmax-1) {  T[j,k] <- 1 - phi(A + d[k] + z[j]) }
  }
}                                     # *** PROGRAM ENDS                                                           

  filename <- "psoriasisREexample.bug"
  ## write model file:
  write.model(psoriasisRE.model,filename)
  ## and let's take a look:
  file.show(filename)


ns<-16 #no of studies 
nt<-8  #no of treatments
Cmax<-4 # max number of categories
meanA=1.096
precA=129
t<-as.matrix(rin[,2:4])
na<-rin[,5]
nc<-rin[,6]
C<-as.matrix(rin[,7:10])
r<-array(c(rin$r11,rin$r12,rin$r13,rin$r21,rin$r22,rin$r23,rin$r31,rin$r32,rin$r33),dim=c(16,3,3))
n<-array(c(rin$n11,rin$n12,rin$n13,rin$n21,rin$n22,rin$n23,rin$n31,rin$n32,rin$n33),3))

data<-list("ns"=ns,"nt"=nt,"Cmax"=Cmax,"meanA"=meanA,"precA"=precA,"t"=t,"na"=na,"nc"=nc,"C"=C,"r"=r,"n"=n)



# Initial values
inits<-list(list(d = c(NA,0),A=1,mu = c(0,sd = 1,z.aux=c(NA,0.66,1.3)),list(d = c(NA,1,2,sd = 1.5,0.5,1)),list(d= c(NA,1),sd = 3,0.1,2)))

PsoriasisRE.sim <- bugs(data,inits,model.file = "C:/Users/User1/Autotask Workplace/Axis Working Folder/Niamh Russell Folder/WinBUGS/psoriasisREexample.bug",parameters = c("d","sd","z","A"),n.chains = 3,n.iter = 50000,debug=TRUE,bugs.directory = "C:/Users/User1/Documents/Winbug/winbugs14_full_patched/winbugs14")

然而，当我尝试运行代码时，它挂在 R 中，并且 Winbugs 不会尝试加载初始化。该过程似乎在编译步骤后停止，并且将消息“使用未定义的节点 n”写入日志。

n，然而，在 R 中被定义为一个 3 维数组，并构成传递给 Winbugs 的数据的一部分。

有人可以帮忙吗。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

r r r2winbugs