dmnorm 错误 "Invalid parent values" -- 从 JAGS 中的 GeoBUGS 手册重现苏格兰唇癌示例

dmnorm error "Invalid parent values" -- Reproducing Scottish Lip Cancer example from the GeoBUGS manual in JAGS

我正在尝试从 JAGS 的 GeoBUGS 手册中复制苏格兰唇癌示例。但是,我不断从 dmnorm 函数 "Invalid parent values" 中收到以下错误。我手动将精度矩阵输入 JAGS,我知道这是可以接受的。不确定发生了什么,如果有新的眼睛看一下就好了。下面是(希望)一个可重现的例子。

#######################################
## Understanding CAR prior with JAGS ##
#######################################

#Load libraries
  library(rjags)
  library(rgdal)
  library(spdep)
  library(utils)

#Download data
  setwd(tempdir())
  setInternet2(use=TRUE)
  download.file('https://geodacenter.org/downloads/data-files/scotlip.zip','scotlip.zip')
  unzip('scotlip.zip')
  scotland = readOGR('scotlip.shp',layer='scotlip')

#Extract data (get rid of islands for simplicity)
  scot.data = scotland@data[which(!1:56%in%c(6,8,11)),c('DISTRICT','CANCER','CEXP','AFF')]

#Extract adjacency matrix (get rid of islands for simplicity)
  scot.adj = poly2nb(scotland[which(!1:56%in%c(6,8,11)),])

#Visualize
  plot(scotland[which(!1:56%in%c(6,8,11)),],border="grey")
  plot(scot.adj,coordinates(scotland[which(!1:56%in%c(6,8,11)),]),pch=16,add=TRUE)

#Construct weight matrix (or proximity matrix)
  scot.weights = nb2mat(scot.adj,style="B")

#Set spatial correlation for ICAR and compute precision matrix
  alpha = 1/max(eigen(scot.weights)$values)
  tau2 = 1
  P = diag(rowSums(scot.weights))%*%(diag(nrow(scot.weights))-alpha*scot.weights)/tau2

#Is P valid?
  library(MASS)
  Sigma = chol2inv(P)
  mvrnorm(n=1,mu=rep(0,dim(Sigma)[1]),Sigma=Sigma) #Useful for simulations later on

#Load the data
  lips.data = list(
    N=nrow(scot.data),
    O=scot.data$CANCER,
    E=scot.data$CEXP,
    X=scot.data$AFF,
    Tau=P
  )

#Model file
  sink("ScottishLipCancer.txt")
  cat("data{

     for(i in 1:N){ zeros[i] <- 0 }

   }model{

      #Likelihood
        for(i in 1:N){
          O[i] ~ dpois(mu[i])
          log(mu[i]) <- log(E[i]) + alpha0 + alpha1 * X[i]/10 + b[i]

          #Area-specific relative risk (for maps)
            RR[i] <- exp(alpha0 + alpha1 * X[i]/10 + b[i])
        }

      #ICAR prior distribution for random effects:
        b[1:N] ~ dmnorm(zeros[],Tau[,])

      #Other priors:
        alpha0 ~ dnorm(0.0,1.0E-5)
        alpha1 ~ dnorm(0.0,1.0E-5)
        tau ~ dgamma(0.5,0.0005)             
        sigma <- sqrt(1/tau)            
        b.mean <- sum(b[])

    }
  ",fill=TRUE)
  sink()

#Initial values for MCMC chains
  lips.inits = list(tau=1,alpha0=0,alpha1=0,b=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))

#Run JAGS
  m = jags.model(file='ScottishLipCancer.txt',data=lips.data,inits=lapply(1:3,function(x)lips.inits),n.chains=3,n.adapt=10)

我想使用 JAGS 的一个挫折是你必须指定联合分布(通过 dmnorm)而不是使用一堆条件,比如 ICAR 之前的 BUGS。任何帮助,将不胜感激!

我个人认为是您的 zeros[] 代码不正确。 jags 无法确定 zeros[] 的大小,因此应该在 'data' 部分中明确说明,或者将其作为 R

中的变量传入

尝试运行

cat("
  data zeros[N];

  model{
    for(i in 1:N){
      zeros[i] <- 0
      }
    ...
  }
")

Martyn Plummer 更正了我的代码中的一些错误并在 JAGS 站点上发布了答案。可以在此处找到该讨论 http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/b862b65b/.