1

我正在尝试学习如何在 WinBUGS 中拟合层次模型。作为练习，我在 R 中模拟了一些具有 3 个级别的剂量响应数据，然后我尝试将其输入 WinBUGS。

模拟设置：

  • 3 级：CENTER, SUBJECT,RESIDUAL
  • 20个中心
  • 每个中心有20个科目
  • 每个受试者都接受了 9 个剂量水平（忽略此设计的现实方面）
  • 9 个剂量水平是：1, 2, 5, 10, 25, 50, 100, 200, 和500
  • 对数（受试者i在剂量下的反应j）=log(exp(l.E0) + exp(l.EMAX)*(Dose/( exp(l.ED50) + Dose))) + residvar
  • l.E0= 0.693（固定效应）
  • l.ED50= l.ED50.MEAN + RANDOM.SUBJECT.ED50 （所以一个主题特定的 l.ED50 效果）
  • l.EMAX= l.EMAX.MEAN + RANDOM.CENTER.EMAX + RANDOM.SUBJECT.EMAX （所以既是中心又是主题特定的 l.EMAX 效果）
  • l.ED50.MEAN= 3.22（固定效应）
  • l.EMAX.MEAN= 5.706（固定效果）
  • RANDOM.CENTER.EMAX ~ N(0, 1.14)
  • RANDOM.SUBJECT.ED50 + RANDOM.SUBJECT.EMAX ~ MVN([0,0], [0.1036, 0.156, 0.156, 0.325]) （所以双变量正态分布，相关性为 0.85，方差分别为 0.1036 和ED500.325 EMAX）

下面我给出了我尝试过的 WinBUGS 代码和数据。该模型在语法上是正确的，但是在编译数据后，我收到了一些我不理解的错误消息。

我想知道是否有人可以告诉我：

  • 如何更正我的代码
  • 或者如果我对多级多元分布进行编码的尝试是错误的，并且是否有更好的方法来做到这一点。

提前感谢您的任何反馈。

克鲁，罗尔

model {
  for (k in 1:nobs) {
    #model for observed data
    yobs.log.prop[k] ~ dnorm(yobs.log.prop.hat[k], tau)          
    # likelihood        
    yobs.log.prop.hat[k] <- log(exp(mu.l.E0) + exp(l.ID.EMAX[nsub.unique.from1[k]]) * (dose[k]/(dose[k] + exp(l.ID.ED50[nsub.unique.from1[k]]))))
    l.ID.EMAX[nsub.unique.from1[k]] <- mu.l.EMAX + u.EMAX.CENTER[ncenter[k]] + u.ID[nsub.unique.from1[k], 1]
    l.ID.ED50[nsub.unique.from1[k]] <- mu.l.ED50 + u.ID[nsub.unique.from1[k], 2]
  }
  for (j in 1:20) {
    #intercenter variability: level1
    u.EMAX.CENTER[j] ~ dnorm(0, tau.l.EMAX)
    u.EMAX.CENTER.pred[j] ~ dnorm(0, tau.l.EMAX)
  }
  for (i in 1:400) {
    #intersubject variability: level 2
    u.ID[i, 1:2] ~ dnorm(zero2[1:2], tau.u2[1:2, 1:2])
  }
  #Priors fixed effects
  mu.l.E0 ~ dunif(-2, 5)
  mu.l.EMAX ~ dunif(1, 10)
  mu.l.ED50 ~ dunif(-2, 5)
  tau <- 1/(sigma*sigma)
  sigma ~ dunif(0,1000)
  #Prior for CENTER effect
  tau.l.EMAX <- 1/(omega.l.EMAX * omega.l.EMAX)
  omega.l.EMAX ~ dunif(0, 100)
  #omega.l.EMAX ~ dnorm(0,0.01)I(0,)
  #Prior for ID effect
  for (m in 1:2) {
    zero2[m] <- 0
  }
  tau.u2[1:2, 1:2] ~ dwish(R[1:2,1:2],2) 
  sigma.u2[1:2, 1:2] <- inverse(tau.u2[1:2, 1:2])
  rho <- sigma.u2[1, 2]/sqrt(sigma.u2[1,1]*sigma.u2[2, 2])
  #omega.l.EMAX~dnorm(0,0.01)I(0,)
}

我想包含数据，但收到了我习惯了很多字符的消息。因此，我只显示开头（带有 Wishart R 矩阵的规范）和我的数据的一个子集。

list(nobs=3600, R = structure(.Data = c(1, 0.85, 0.85, 1), .Dim = c(2, 2)),
     yobs.log.prop=c(4.29774919346511, ...), #... are the remaining 3599 responses
     dose=c(1,  ...), #... are the remaining 3599 dose values
     ncenter=c(1, ....), #... are the remaining 3599 center identifications
     nsub.unique.from1=c(1...)) #... are the remaining 3599 subject identifications
4

0 回答 0