1

我正在尝试复制 JAGS 中 GeoBUGS 手册中的苏格兰唇癌示例。但是，我不断从 dmnorm 函数“无效的父值”中收到以下错误。我手动将精度矩阵输入 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 这样的条件。任何帮助，将不胜感激！

4

2 回答 2

0

就个人而言，我认为这是您的 zeros[] 代码不正确。jags 无法判断 zeros[] 的大小，因此应该在“数据”部分中明确说明，或者将其作为 R 中的变量传递

尝试运行

cat("
  data zeros[N];

  model{
    for(i in 1:N){
      zeros[i] <- 0
      }
    ...
  }
")
于 2015-02-22T15:37:39.533 回答
0

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

于 2015-03-02T21:40:40.163 回答