我无法让 Jags 中的“排名”功能正常工作。下面是改编自 Winbugs 的代码、模型和数据。不起作用的一点是：
myorder[i] <- rank(aux.u[],i)
，这会引发以下错误：
RUNTIME ERROR:
Incorrect number of parameters in function rank
我知道 Jags 中的功能等级不同。所以，当我尝试使用 somwhing 重新编码它时：
myorderx <- rank(aux.u[]);myorder[i]<-myorderx[i];
我得到错误：
Attempt to redefine node myorderx[1:4]
任何建议表示赞赏。非常感谢，丽索
模型
rm(list = ls())
library(coda)
library(rjags)
library(pracma);
library(nmathresh)
setwd("/Users/test")
nummodel <- function(){
for(i in 1:NS) {
w[i,1] <- 0
delta[i,t[i,1]] <- 0
mu[i] ~ dnorm(0,.0001) # Fixed study effects
for (k in 1:na[i]) {
r[i,k] ~ dbin(p[i,t[i,k]],n[i,k]) # Binomial likelihood for data
logit(p[i,t[i,k]])<-mu[i] + delta[i,t[i,k]] } # Model for log odds parameters
for (k in 2:na[i]) {
delta[i,t[i,k]] ~ dnorm(md[i,t[i,k]],precd[i,t[i,k]]) # Random treatment effects
md[i,t[i,k]] <- d[t[i,k]] - d[t[i,1]] + sw[i,k] # Mean of LOR distributions
precd[i,t[i,k]] <- prec[t[i,1],t[i,k]]*2*(k-1)/k # Precision of LOR distributions
w[i,k] <- (delta[i,t[i,k]] - d[t[i,k]] + d[t[i,1]]) # Adjustment, multi-arm RCTs
sw[i,k] <-sum(w[i,1:(k-1)])/(k-1) } # Cumulative adjustment, multi-arm RCTs
}
d[1]<-0
for (k in 2:NT) {d[k] ~ dnorm(0,.0001) } # Vague priors for basic parameters
for(j in 1:3){
prec[j,j]<-1
for(k in (j+1):4){
prec[j,k]<-1/tausq[j,k]
prec[k,j]<-prec[j,k] }}
for(k in 1:4){
v.a[k]~dlnorm(-3.31,0.346) # Informative prior chosen to imply (approximately) chosen data-based prior for contrast heterogeneity variance
sd.a[k]<-sqrt(v.a[k])
}
pi.half<-1.5708
for(i in 1:3) {for(j in (i+1):4){
g[j,i]<-0
tausq[i,j]<-v.a[i]+v.a[j]-2*rho.star[i,j]*sd.a[i]*sd.a[j]
tau[i,j] <- sqrt(tausq[i,j]) }}
# Implementing random permutation
for(i in 1:4) { for(j in 1:4) {
rho[i,j]<-inprod(g[, i],g[, j])
rho.star[i,j]<- rho[ myorder[i], myorder[j] ] }}
for(i in 1:4) {
aux.u[i] ~ dunif(0, 1)
myorder[i] <- rank(aux.u[],i)
#myorderx <- rank(aux.u[]);myorder[i]<-myorderx[i];
}
# Constructing entries of upper-triangular matrix for Cholesky decomposition
g[1,1]<-1
g[2,2]<-sin.a[1,2]
g[3,3]<-sin.a[1,3]*sin.a[2,3]
g[4,4]<-sin.a[1,4]*sin.a[2,4]*sin.a[3,4]
g[1,2]<-cos.a[1,2]
g[1,3]<-cos.a[1,3]
g[1,4]<-cos.a[1,4]
g[2,3]<-sin.a[1,3]*cos.a[2,3]
g[2,4]<-sin.a[1,4]*cos.a[2,4]
g[3,4]<-sin.a[1,4]*sin.a[2,4]*cos.a[3,4]
# Beta prior for cos(a[i,j])
for (i in 1:3) {
for (j in (i+1):4) {
cos.a[i,j] ~ dbeta(0.93,1.07) # From Table 3
sin2.a[i,j] <- 1 - pow(cos.a[i,j] , 2)
sin.a[i,j] <- pow(sin2.a[i,j] , 1/2) } }
# Defining pairwise ORs
for (c in 1:(NT-1))
{ for (k in (c+1):NT)
{ lor[c,k] <- d[k] - d[c]
log(or[c,k]) <- lor[c,k] }
}
}
write.model(nummodel, "Model.txt")
model.file1 = paste(getwd(),"Model.txt", sep="/")
为 JAGS 定义一些 MCMC 参数
nchains <- 3; # How Many Chains?
nadapt<-200
nburnin <- 200; # How Many Burn-in Samples?
nsamples <- 300; # How Many Recorded Samples?
nthin <- 20;
数据
r1<-c(9,11,75,2,58,0,3,1,6,79,18,64,5,20,0,8,95,15,78,69,20,7,12,9)
n1<-c(140,78,731,106,549,33,100,31,39,702,671,642,62,234,20,116,1107,187,584,1177,49,66,76,55)
r2<-c(23,12,363,9,237,9,31,26,17,77,21,107,8,34,9,19,143,36,73,54,16,32,20,3)
n2<-c(140,85,714,205,1561,48,98,95,77,694,535,761,90,237,20,149,1031,504,675,888,43,127,74,26)
r3<-c(10,29,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
n3<-c(138,170,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
t1<-c(1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,3,3)
t2<-c(3,3,3,3,3,3,3,3,3,2,2,3,3,3,4,2,3,3,3,3,3,4,4,4)
t3<-c(4,4,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
na<-c(3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2)
datasrc <- data.frame(r1,r2,r3,n1,n2,n3,t1,t2,t3,na)
datastruct<-list(NS=length(datasrc$t1),NT=18,na=c(datasrc$na),
t=structure(.Data=c(datasrc$t1,datasrc$t2,datasrc$t3), .Dim=c(length(datasrc$t1),3)),
r=structure(.Data=c(datasrc$r1,datasrc$r2,datasrc$r3), .Dim=c(length(datasrc$t1),3)),
n=structure(.Data=c(datasrc$n1,datasrc$n2,datasrc$n3), .Dim=c(length(datasrc$t1),3)))
运行模型
parameters = c("lor");
model.file1 = "/Users/test/Model.txt"
mod1 <- jags.model(file =model.file1, data=datastruct, n.chains=nchains, n.adapt=nadapt);