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/.
我正在尝试从 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/.