How to make leave-one-out cross-validation using R2jags or rjags

128 views Asked by At

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 :)

0

There are 0 answers