我在使用以下代码编译 Jags 模型时遇到问题。似乎正在尝试重新定义参数 delta （第 17 行： delta[i, si[i, k], hh] ~ dnorm(md[i, si[i, k], hh], taud[i, si[i, k], hh])) 在模型内第三个循环的开头，但我不确定为什么？
型号为：
model
{
for (hh in 1:nrepeat) {
for (i in 1:ns) {
w[i, 1, hh] <- 0.00000E+00
j[i, 1, hh] <- 0.00000E+00
delta[i, bi[i, hh], hh] <- 0.00000E+00
mu[i, hh] ~ dnorm(0.00000E+00, 1.00000E-05)
beta0[i, hh] ~ dnorm(0.00000E+00, 1.00000E-05)
b1[i, 1, hh] <- 0.00000E+00
for (k in 1:na[i]) {
index[i, k, hh] <- split[i, hh] * (equals(t[i,
k], pair[hh, 1]) + equals(t[i, k], pair[hh,
2]))
}
for (k in 2:na[i]) {
delta[i, si[i, k], hh] ~ dnorm(md[i, si[i, k],
hh], taud[i, si[i, k], hh])
md[i, si[i, k], hh] <- (d[si[i, k], hh] - d[bi[i,
hh], hh] + sw[i, k, hh] + b1[i, k, hh]) * (1 -
index[i, m[i, k], hh]) + direct[hh] * index[i,
m[i, k], hh]
j[i, k, hh] <- k - (equals(1, split[i, hh]) *
step(k - 3))
taud[i, si[i, k], hh] <- tau[hh] * 2 * (j[i,
k, hh] - 1)/j[i, k, hh]
w[i, k, hh] <- (delta[i, si[i, k], hh] - ((d[si[i,
k], hh] - d[bi[i, hh], hh]) + b1[i, k, hh])) *
(1 - index[i, k, hh])
sw[i, k, hh] <- sum(w[i, 1:(k - 1), hh])/(j[i,
k, hh] - 1)
b1[i, k, hh] <- (beta[t[i, k], hh] * (x[i] -
mx))
}
}
for (l in 1:np) {
r[l] ~ dbin(p[l], n[l])
logit(p[l]) <- mu[s[l], hh] + beta0[s[l], hh] * (x[l] -
mx) + delta[s[l], tipd[l], hh] + (deltab[l, hh] *
(1 - equals(tipd[l], bi[s[l], hh])))
rhat[l] <- p[l] * n[l]
dev[l] <- 2 * (r[l] * (log(r[l]) - log(rhat[l])) +
(n[l] - r[l]) * (log(n[l] - r[l]) - log(n[l] -
rhat[l])))
index2[l, hh] <- split[s[l], hh] * (equals(tipd[l],
pair[hh, 1]) + equals(tipd[l], pair[hh, 2]))
deltab[l, hh] <- (beta[tipd[l], hh] - beta[bi[s[l],
hh], hh]) * (x[l] - mx) * (1 - index2[l, hh]) +
directbeta[hh] * (x[l] - mx) * (index2[l, hh])
}
totresdev[hh] <- sum(dev[])
direct[hh] ~ dnorm(0.00000E+00, 1.00000E-05)
directbeta[hh] ~ dnorm(0.00000E+00, 1.00000E-05)
d[1, hh] <- 0.00000E+00
beta[1, hh] <- 0.00000E+00
sd[hh] ~ dunif(0.00000E+00, 10)
tau[hh] <- pow(sd[hh], -2)
tausq[hh] <- sd[hh] * sd[hh]
for (k in 2:nt) {
d[k, hh] ~ dnorm(0.00000E+00, 1.00000E-05)
beta[k, hh] ~ dnorm(0.00000E+00, 1.00000E-05)
}
for (k in 1:nt) {
for (v in 1:nz) {
dz[v, k, hh] <- d[k, hh] - (beta[k, hh]) * (mx -
z[v])
}
}
for (c in 1:(nt - 1)) {
for (k in (c + 1):nt) {
betas[c, k, hh] <- beta[k, hh] - beta[c, hh]
lor[c, k, hh] <- (d[k, hh] - d[c, hh])
for (v in 1:nz) {
lorz[v, c, k] <- (dz[v, k, hh] - dz[v, c, hh])
}
}
}
for (v in 1:nz) {
directz[v, hh] <- direct[hh] - (directbeta[hh]) *
(mx - z[v])
directorz[v, hh] <- exp(directz[v, hh])
}
for (v in 1:nz) {
diff[v, hh] <- directz[v, hh] - lorz[v, pair[hh,
1], pair[hh, 2]]
prob[v, hh] <- step(diff[v, hh])
}
}
}
数据文件可以在这里获取： 数据和模型文件
JAGS 的输入：
rm(list = ls())
library(coda)
library(rjags)
library(netmeta);
#LOAD FUNCTIONS TO SHAPE DATA
#CHECK IF PAIR(X,Y) IN ROW I OF DATA AND GIVE BASELINE FOR DATA ROW I
PairXY <- function(treat, pair)
{
N <- nrow(treat)
out <- cbind(split=rep(0,N), b=rep(0,N))
for (i in 1:N) {
pos <- match(pair, treat[i,], nomatch=0) # lenght = length(pair) = 2
out[i,1] <- ifelse(prod(pos)>0, 1, 0) # 1 if pair in line i, 0 o.w.
out[i,2] <- ifelse(prod(pos)==0, 1, pos[1])
}
out
}
# GIVES NA-1 INDEXES TO SWEEP NON-BASELINE ARMS ONLY
NonbaseSweep <- function(index, na)
{
N <- NROW(na)
C <- max(na)
out <- matrix(nrow=N, ncol=C)
for (i in 1:N) {
for (k in 2:na[i]) {
out[i,k] <- k - (index[i,"b"] >= k)
}
}
out
}
# BUILDS MATRIX WITH NON-BASELINE TREATMENTS
Sweeptreat <- function(treat, m)
{
N <- NROW(treat)
C <- NCOL(m)
out <- matrix(nrow=N, ncol=C)
for (i in 1:N) {
for (k in 2:C) {
out[i,k] <- treat[i,m[i,k]]
}
}
out
}
## BUILDS VECTOR WITH BASELINE TREATMENTS
Basetreat <- function(treat, b)
{
N <- nrow(treat)
out <- rep(0,N)
for (i in 1:N) {
out[i] <- treat[i,b[i]]
}
out
}
# Data
datasrc <- read.csv('OR_Responders.csv', header=TRUE)
datasrc[] <- lapply(datasrc, function(x){
x[is.nan(x)] <- NA
x
})
datasrc<-subset(datasrc,datasrc$XSev1!="NA")
mx<-mean(c(datasrc$XSev1,datasrc$XSev2,datasrc$XSev3),na.rm = TRUE);
dat1<-datasrc
dat1$study<-1:length(dat1$t1)
dat2 <- pairwise(list(t1, t2, t3),event=list(r1,r2,r3),n=list(n1, n2, n3),data=dat1, studlab=study)
na=dat1$na #NUMBER OF ARMS IN EACH STUDY
t=cbind(dat1$t1,dat1$t2,dat1$t3, deparse.level = 0) #TREATMENT NUMBER
s=dat2$studlab #STUDY NUMBER
#y=dat2$TE #OUTCOME
r=dat2$event1+dat2$event2
n=dat2$n1+dat2$n2
tipd=dat2$treat2 #TREATMENT NUMBER
x=dat2$XSevag #COVARIATE VALUES
#x=dat2$age/12 #COVARIATE VALUES
ns=max(s) #NMUBER OF TRIALS
nt=22 #NUMBER OF TREATMENTS
np=length(x) #NUMBER OF PATIENTS
mx=mean(x) #AVERAGE COVARIATE VALUE
z=c(10,20,30,40,50) #CHOSEN COVARIATE VALUES AT WHICH TREATMENT EFFECTS ARE REQUIRED TO BE ESTIMATED.
nz=length(z) #NUMBER OF CHOSEN COVARIATE VALUES
#CHOOSE NODE TO SPLIT
parameters = c("prob","diff");model.file1 = "ModelC1.txt"
contre <- read.csv('Treat.csv', header=TRUE)
pair<-contre[,2:3];pair<-data.matrix(pair, rownames.force = NA)
split <- array(0, dim=c(ns,length(contre[,1])))
bi <- array(0, dim=c(ns,length(contre[,1])))
for (i in 1:length(contre[,1])) {
pairx <- c(contre[i,2],contre[i,3]);
checkPair <- PairXY(t, pairx);split[,i]<-checkPair[,"split"]
bix <- Basetreat(t, checkPair[,"b"]);bi[,i] = bix;
m <- NonbaseSweep(checkPair, na);#m[,i] = mtt;
si <- Sweeptreat(t,m);#si[,i] = six;
}
# Defining some MCMC parameters for JAGS
nchains <- 3; # How Many Chains?
nadapt<-100
nburnin <- 100; # How Many Burn-in Samples?
nsamples <- 100; # How Many Recorded Samples?
nthin <- 10;
datastruct<-list(z=z, nz=nz,ns=ns,nt=nt,mx=mx,s=s,na=na,nrepeat=length(pair[,1]),split = split,m=m,bi=bi,si=si,pair=pair,x=x,tipd=tipd,t=t,r=r,n=n,np=np)
mod1 <- jags.model(file =model.file1, data=datastruct, n.chains=nchains, n.adapt=nadapt);
# update(mod1, nburnin);
# sim<-coda.samples(mod1, parameters, nsamples, thin = nthin, na.rm=TRUE);
# simsum<-summary(sim);
错误信息是：
Compiling model graph
Resolving undeclared variables
Allocating nodes
Deleting model
Error in jags.model(file = model.file1, data = datastruct, n.chains = nchains, :
RUNTIME ERROR:
Compilation error on line 17.
Attempt to redefine node delta[15,1,1]
任何提示将不胜感激。谢谢丽索