У меня проблема с компиляцией модели Jags, преобразованной из модели Openbugs, где она работает нормально, с использованием следующего кода.
Кажется, что параметр w нельзя индексировать в конце второго цикла в модели, но я не уверен, почему?
rm(list = ls())
library(rjags)
library(jagsUI)
My=getwd()
cat('
model{ # *** PROGRAM STARTS
for(i in 1:ns){ # LOOP THROUGH STUDIES
w[i,1] <- 0 # adjustment for multi-arm trials is zero for control arm
delta[i,1] <- 0 # treatment effect is zero for control arm
mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines
for (k in 1:na[i]) { # LOOP THROUGH ARMS
r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood
# model for linear predictor, covariate effect relative to treat in arm 1 (centring)
logit(p[i,k]) <- mu[i] + delta[i,k] + (beta[t[i,k]]-beta[t[i,1]]) * (x[i]-mx)
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators
dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k])))
}
resdev[i] <- sum(dev[i,1:na[i]]) # summed residual deviance contribution for this trial
for (k in 2:na[i]) { # LOOP THROUGH ARMS
delta[i,k] ~ dnorm(md[i,k],taud[i,k]) # trial-specific LOR distributions
md[i,k] <- d[t[i,k]] - d[t[i,1]] + sw[i,k] # mean of LOR distributions (with multi-arm trial correction)
taud[i,k] <- tau *2*(k-1)/k # precision of LOR distributions (with multi-arm trial correction)
w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]]) # adjustment for multi-arm RCTs
sw[i,k] <- sum(w[i,1:k-1])/(k-1) # cumulative adjustment for multi-arm trials
}
}
totresdev <- sum(resdev[]) # Total Residual Deviance
d[1]<-0 # treatment effect is zero for reference treatment
beta[1] <- 0 # covariate effect is zero for reference treatment
for (k in 2:nt){ # LOOP THROUGH TREATMENTS
d[k] ~ dnorm(0,.0001) # vague priors for treatment effects
beta[k] <- B # common covariate effect
}
B ~ dnorm(0,.0001) # vague prior for covariate effect
sd ~ dunif(0,5) # vague prior for between-trial SD
tau <- pow(sd,-2) # between-trial precision = (1/between-trial variance)
} # *** PROGRAM ENDS
', file={f <- tempfile()})
data.src <- read.table(textConnection('
t1 t2 na n1 n2 r1 r2 x id #Study name
1 3 2 63 65 9 28 6.85 1 #Kim 2007 (37)
1 3 2 200 207 19 81 10.95 2 #DE019 Trial (36)
1 3 2 62 67 5 37 11.65 3 #ARMADA Trial (34)
1 2 2 199 393 15 146 6.15 4 #RAPID 1 Trial (40)
1 2 2 127 246 4 80 5.85 5 #RAPID 2 Trial (41)
1 5 2 363 360 33 110 8.10 6 #START Study (57)
1 5 2 110 165 22 61 7.85 7 #ATTEST Trial (51)
1 5 2 47 49 0 15 8.30 8 #Abe 2006 (50)
1 4 2 30 59 1 23 13.00 9 #Weinblatt 1999 (49)
1 6 2 40 40 5 5 11.25 11 #Strand 2006 (62)
1 7 2 49 50 14 26 0.92 12 #CHARISMA Study (64)
1 7 2 204 205 22 90 7.65 13 #OPTION Trial (67)
'), header=TRUE)
data<-list(ns=12, nt=7,mx=8.21,na=c(data.src$na),x=c(data.src$x),
t=structure(.Data=c(data.src$t1,data.src$t2), .Dim=c(12,2)),
r=structure(.Data=c(data.src$r1,data.src$r2), .Dim=c(12,2)),
n=structure(.Data=c(data.src$n1,data.src$n2), .Dim=c(12,2)))
inits <- function(){
#chain 1
list(d=c( NA, 0,0,0,0,0,0), mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0), sd=1, B=0)
#chain 2
list(d=c( NA, -1,1,-1,1,-1,1), mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3), sd=0.5, B=-1)
#chain 3
list(d=c( NA, 2,-2,2,-2,2,-2), mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0, 5, 0), sd=3, B=5)
}
parameters = c("beta","d")
sim <- jags(data = data,inits = inits,parameters.to.save = parameters,model.file = f,n.chains = 3,n.adapt = 100,n.iter = 1000,n.burnin = 500,n.thin = 2)
Я получаю следующее сообщение об ошибке:
Processing function input.......
Done.
Compiling model graph
Resolving undeclared variables
Allocating nodes
Deleting model
Error in jags.model(file = model.file, data = data, inits = inits, n.chains = n.chains, :
RUNTIME ERROR:
Compilation error on line 22.
Index out of range taking subset of w
Любая подсказка будет оценена.
Спасибо Лисо