I'm trying to make leave-one-out cross-validation using R2jags. I do it easily when I use R2WinBUGS instead with the same codes, but jags() give me this error :
Compiling model graph
Resolving undeclared variables
Allocating nodes
Deleting model
Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
RUNTIME ERROR:
Compilation error on line 18.
Y[1,1:2] is partly observed and partly missing
Here are the codes and model I use
model
model_DH <- "model{
se[index,2]~dunif(0.0001,15) #prior on the missing SE in the validation study i=index
pred.delta2<-delta[index,2] #predicted effect on the final outcome in the validation study
#within study precision matrix
#rho_w<-corr_w #if known, a fixed correlation can be used
rho_w~dunif(0,0.999) #we assume positive within-study correlation
for (i in 1:num) {
Prec_w[i,1:2,1:2] <- inverse(sigma[i,1:2,1:2]) # variance-covariance matrix
#covariance matrix for the j-th study
sigma[i,1,1]<-pow(se[i,1],2)
sigma[i,2,2]<-pow(se[i,2],2)
sigma[i,1,2]<-sqrt(sigma[i,1,1])*sqrt(sigma[i,2,2])*rho_w
sigma[i,2,1]<-sqrt(sigma[i,1,1])*sqrt(sigma[i,2,2])*rho_w
}
# Daniels and Hughes formulation for the between-studies model:
for (i in 1:num) {
Y[i,1:2]~dmnorm(delta[i,1:2], Prec_w[i,1:2,1:2])
# product normal formulation for the between study part:
delta[i,1]~dnorm(0.0, 0.001)
delta[i,2]~dnorm(eta2[i],prec2)
eta2[i]<-lambda0+lambda1*delta[i,1]
}
lambda0~dnorm(0.0, 1.0E-3)
lambda1~dnorm(0.0, 1.0E-3)
psi.2~dunif(0,2)
psi2.sq<-pow(psi.2,2)
prec2<-1/psi2.sq
}"
writeLines(model_DH,"model_name.txt")
library(R2jags)
data preparation
num <- 20
trials <- as.character(letters[1:20])
set.seed(123)
y1 <- rnorm(num)
y2 <- rnorm(num, mean=2, sd=0.56)
se1 <- seq(0.1, 0.6, length.out=num)
se2 <- seq(0.2, 0.7, length.out=num)
# parameters
n.iter<- 100000 # Number of iterations
n.burnin<-50000 # Burn-in
pred.delta2.DaHu <- pred.cri.DaHu <- bias.DaHu <- cil.DaHu <- ciu.DaHu <- sd.DaHu <-array(0,num)
for(i in 1:num){
index<-i
Y<- se <- array(0,dim=c(num,2))
Y[,1] <- y1
Y[,2] <- y2
se[,1] <- se1
se[,2]<-se2
#setting final outcome in the validation study i to misisng value to be predicted
Y[i,2]<-NA
se[i,2]<-NA
data=list(Y=Y, se=se, num=num, index=i)
#inits #setting initial values for missing final outcome in the validation study
Y0 <- se0 <-structure(rep(NA,num*2),dim=c(num,2))
Y0[index,2] <- -0.5
se0[index,2] <- 0.2
inits <- list(list(psi.2 = 0.25, lambda1 = 0.0, lambda0 = 0.0, rho_w=0.25,
delta = structure(.Data = c(rep(-0.5,num),rep(-0.5,num)),
.Dim=c(num,2)), Y=Y0,se=se0))
#monitoring
para=c("psi2.sq", "pred.delta2")
fit.Dahu<-jags(data=data,
inits=inits, parameters.to.save=para,
model.file="model_name.txt",
n.chains=1, n.iter=n.iter, n.burnin=n.burnin,
DIC=F, working.directory=getwd(),
jags.seed = 123, digits = 3)
x <- fit.DaHu$BUGSoutput$sims.matrix[,grep("pred.delta2",colnames(fit.DaHu$BUGSoutput$sims.matrix))]
pred.delta2.DaHu[i]<-mean(x)
bias.DaHu[i]<-logHR_OS[i]-pred.delta2.DaHu[i]
sd.DaHu[i] <-fit.DaHu$BUGSoutput$summary[grep("pred.delta2",rownames(fit.DaHu$BUGSoutput$summary)),c(2)]
pred.cri.DaHu[i]<-sqrt(logHR_OS_se[i]^2+sd.DaHu[i]^2)*2.00*1.96
cil.DaHu[i] <- pred.delta2.DaHu[i]-pred.cri.DaHu[i]/2.0
ciu.DaHu[i] <- pred.delta2.DaHu[i]+pred.cri.DaHu[i]/2.0
cat("\n End of ", i, "eme iteration \n")
}
Thanks a lot for your answers :)