doseresNMA:

#!!!!!!! ask Georgia:   random/fixed effect to which should be applied: delta, beta1, beta2, gamma, class.b
#!!!!!!! ask Georgia: in case of metareg, add the option to either indepndent (bb[bb[k] ~ dnorm(0,0.001)]) or common covariate effect ( bb[k] ~ dnorm(b,precbb))
#!!!!!!! I added the absolute effect part but make it nicer
#!!!!!!! organize it better
# check Gerogia, that is not true? isn't it?
# for (i in 1:ndrugs) {
#   for (j in 1:ndrugs) {
#     bb1[i,j] <- b1[i]-b1[j]
#     bb2[i,j] <- b2[i]-b2[j]
#   }
# }

makeJAGSmodel <- function(doselink,metareg,beta.effects,consistency,class.effect){
  if(class.effect){
    if(doselink=="linear"){
      dose.str <- "d[i,k] <-  (beta1[t[i,k],class[i,k]]*dose1[i,k])-(beta1[t[i,1],class[i,1]]*dose1[i,1])"
      prior.b1.str <- "
      beta1[ref.index,refclass.index] <- 0
       b1[refclass.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
      beta1[k,c] ~ dnorm(b1[c],precb)
    }
  }
   for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
  b1[c]~dnorm(0,0.001)
   }
# Absolute effect
  for (i in 1:ns) { ## for each study
    rr[i,1] ~ dbinom(p0[i],nn[i,1])
    logit(p0[i]) <- z[i]
    z[i] ~ dnorm(Z, prec.z)
  }
  # all comparisons
for (i in 1:ndrugs) {
  for (j in 1:ndrugs) {
    bb1[i,j] <- b1[i]-b1[j]
  }
}
  # priors
  Z ~ dnorm(0, 0.001)
  prec.z <- 1/v.z
  v.z <- sigma.z * sigma.z
  sigma.z ~ dunif(0,10) #!! sprintf not accept dollar sign

  for (k in c(1:(refclass.index-1),(refclass.index+1):nc)){
  for( j in 1:nd.new){
    OR[j,k] <- exp(b1[k]*new.dose[j])
    odds.drug[j,k] <- OR[j,k]*exp(Z)
    p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
  }
  }

"
      prior.b2.str <- ""
    }

    if(doselink== "spline"){
      dose.str <- "d[i,k] <-  (beta1[class[i,k],t[i,k]]*dose1[i,k])-(beta1[class[i,1],t[i,1]]*dose1[i,k]) + (beta2[class[i,k],t[i,k]]*dose2[i,k])-(beta2[class[i,1],t[i,1]]*dose2[i,k])"
      prior.b1.str <- "
      beta1[refclass.index,ref.index] <- 0
      b1[refclass.index] <- 0
      for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
      for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
      beta1[c,k] ~ dnorm(b1[c],precb)
    }
  }

  for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
  b1[c]~dnorm(0,0.001)
  }
   taub~dunif(0,5)      # common standard deviation is given a vague prior
precb<-1/pow(taub,2) # common precision in RE-NMA model
"
      prior.b2.str <- "
      beta2[refclass.index,ref.index] <- 0
      b2[refclass.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
      beta2[c,k] ~ dnorm(b2[c],precb)
    }
  }
  for (c in c(1:(refclass.index-1),(refclass.index+1):nc)) {
  b2[c]~dnorm(0,0.001)
   }

# all comparisons
for (i in 1:nc) {
  for (j in 1:nc) {
    bb1[i,j] <- b1[i]-b1[j]
    bb2[i,j] <- b2[i]-b2[j]
  }
}
   # Absolute effect
  for (i in 1:ns) { ## for each study
    rr[i,1] ~ dbinom(p0[i],nn[i,1])
    logit(p0[i]) <- z[i]
    z[i] ~ dnorm(Z, prec.z)
  }
  # priors
  Z ~ dnorm(0, 0.001)
  prec.z <- 1/v.z
  v.z <- sigma.z * sigma.z
  sigma.z ~ dunif(0,10) #!! sprintf not accept dollar sign

  for (k in c(1:(refclass.index-1),(refclass.index+1):nc)){
  for( j in 1:nd.new){
    OR[j,k] <- exp(b1[k]*new.dose[j]+ b2[k]*f.new.dose[j])
    odds.drug[j,k] <- OR[j,k]*exp(Z)
    p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
  }
  }
"
    }
  }else{
    if(consistency){
    if(doselink=="linear"){
      dose.str <- "d[i,k] <-  (beta1[t[i,k]]*dose1[i,k])-(beta1[t[i,1]]*dose1[i,1])"
      prior.b1.str <- "beta1[ref.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    beta1[k] ~ dnorm(0,0.001)
  }
      for (i in 1:ndrugs) {
  for (j in 1:ndrugs) {
    bb1[i,j] <- beta1[i]-beta1[j]
  }
}"
      prior.b2.str <- ""
    }
    if(doselink== "spline"){
      dose.str <- "d[i,k] <-  ((beta1[i,t[i,k]]*dose1[i,k]) + (beta2[i,t[i,k]]*dose2[i,k]))-((beta1[i,t[i,1]]*dose1[i,1])+(beta2[i,t[i,1]]*dose2[i,1]))"
  #     prior.b1.str <- "
  #     b1[ref.index] <- 0
  # for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
  #   b1[k] ~ dnorm(0,0.001)
  # }"
  #     prior.b2.str <- "
  #           b2[ref.index] <- 0
  # for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
  #   b2[k] ~ dnorm(0,0.001)
  # }

  # all comparisons
# for (i in 1:ndrugs) {
#   for (j in 1:ndrugs) {
#     bb1[i,j] <- b1[i]-b1[j]
#     bb2[i,j] <- b2[i]-b2[j]
#   }
# }
      # Absolute effect
      abs.eff <- " for (i in 1:ns) { ## for each study
    rr[i,1] ~ dbinom(p0[i],nn[i,1])
    logit(p0[i]) <- z[i]
    z[i] ~ dnorm(Z, prec.z)
  }
  # priors
  Z ~ dnorm(0, 0.001)
  prec.z <- 1/v.z
  v.z <- sigma.z * sigma.z
  sigma.z ~ dunif(0,10) #!! sprintf not accept dollar sign

  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
  for( j in 1:nd.new[k]){
    OR[j,k] <- exp(b1[k]*new.dose[k,j]+  b2[k]*f.new.dose[k,j])
    odds.drug[j,k] <- OR[j,k]*exp(Z)
    p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
  }
  }"
    }
    }else{
      if(doselink=="linear"){
        dose.str <- "d[i,k] <-  beta1[t[i,1],t[i,k]]*dose1[i,k]"
        prior.b1.str <- "

  for (i in 1:ncomp){
  beta1[t1[i],t2[i]] <- bb1[t1[i],t2[i]]
      bb1[t1[i],t2[i]] ~ dnorm(0,0.001)
  }"
        prior.b2.str <- ""
      }

      if(doselink== "spline"){
        dose.str <- "d[i,k] <-  beta1[t[i,1],t[i,k]]*dose1[i,k] + beta2[t[i,1],t[i,k]]*dose2[i,k]"
        prior.b1.str <- "
  for (i in 1:ncomp){
  beta1[t1[i],t2[i]] <- bb1[t1[i],t2[i]]
      bb1[t1[i],t2[i]] ~ dnorm(0,0.001)
  }"
        prior.b2.str <- "
  for (i in 1:ncomp){
  beta2[t1[i],t2[i]] <- bb2[t1[i],t2[i]]
      bb2[t1[i],t2[i]] ~ dnorm(0,0.001)
  }"
    }
  }
}
  if (metareg) {
    metareg.str <- "+ (gamma[t[i,k]] - gamma[t[i,1]])*cov[i]"
    prior.b.str <- "
    gamma[ref.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    gamma[k] ~ dnorm(g,prec.g)

  }
    g ~dnorm(0,0.001)
tau.g~dunif(0,5)      # common standard deviation is given a vague prior
prec.g<-1/pow(tau.g,2) # common precision in RE-NMA model
    "

    } else {metareg.str <- ""
  prior.b.str<- ""
  }

if(class.effect){
  beta1.effect <- ""
  beta2.effect <- ""
}else{
  if (beta.effects == "fixed"){
    beta1.effect <- "

     for(i in 1:ns){
     beta1[i,ref.index] <- 0
   for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    beta1[i,k] <- b1[k]
 }}
    "
    beta2.effect <- "
     for(i in 1:ns){
     beta2[i,ref.index] <- 0
   for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
     beta2[i,k] <- b2[k]
   }
    }
    "
  }


  if (beta.effects == "random"){
    beta1.effect <- "
    for(i in 1:ns){
    beta1[i,ref.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    beta1[i,k] ~ dnorm( b1[k],prec.beta)
  }
    }
    tau.beta~dunif(0,5)      # common standard deviation is given a vague prio
    prec.beta<-1/pow(tau.beta,2) # common precision in RE-NMA model
    "
    beta2.effect <- "
    for(i in 1:ns){
    beta2[i,ref.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    beta2[i,k] ~ dnorm(b2[k],prec.beta)
  }
    }"
  }
}
  link.str <- "logit(p[i,k]) <- u[i] + delta[i,k]"
  link.str <- paste0(link.str, metareg.str)

  if(consistency){

    code.str <- sprintf("
    #model{

  for(i in 1:ns){ # Run Through all ns trials
    w[i,1] <-0
    delta[i,1]<-0
    d[i,1] <- 0  #  0 for baseline arms
    u[i] ~ dnorm(0,0.001)

    for (k in 1:na[i]){
      r[i,k] ~ dbin(p[i,k],n[i,k])
      %s
      # deviance
      rhat[i,k] <- p[i,k] * n[i,k]
      dev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (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]])
    for (k in 2:na[i]){
     %s # dose-response fun
      # adjust for within multi-arm correlations
      delta[i,k] <- d[i,k]
      # delta[i,k] ~ dnorm(md[i,k],prec[i,k])
      #
      # md[i,k] <- d[i,k] + sw[i,k]
      # prec[i,k] <- precd *2*(k-1)/k
      # w[i,k] <- delta[i,k] - d[i,k]
      # sw[i,k] <-sum(w[i,1:(k-1)])/(k-1)
    }
  }
%s # beta1 is RE or FE
%s # beta2 is RE or FE

b1[ref.index] <- 0.00000E+00
    for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)) {
      b1[k] ~ dnorm(0.00000E+00, 0.001)
      # b[1:2,k]~dmnorm(b.prior[1:2],inv.R[1:2,1:2])
      # b1[k] <- b[1,k]
      # b2[k] <- b[2,k]
    }
     b2[ref.index] <- 0.00000E+00
    for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)) {
      b2[k] ~ dnorm(0.00000E+00, 0.001)
    }

    # b.prior[1] <- 0
    # b.prior[2] <- 0
    #
    # inv.R ~ dwish(Omega[,], 2)
    # Omega[1,1] <- 1
    # Omega[2,2] <- 1
    #
    # Omega[1,2] <- 0
    # Omega[2,1] <- 0
%s # prior for gamma
%s # absolute effect
# heterogeniety across trials

tau~dunif(0,1)      # common standard deviation is given a vague prior
precd<-pow(tau,-2) # common precision in RE-NMA model
totresdev <- sum(resdev[])
#}
  ",link.str,dose.str,beta1.effect,beta2.effect, prior.b.str,abs.eff)


    }else{
    code.str <- sprintf("

 for(i in 1:ns){ # Run Through all ns trials

    w[i,1] <-0
    delta[i,1]<-0
    u[i] ~ dnorm(0,0.001)

    for (k in 1:na[i]){
      r[i,k] ~ dbin(p[i,k],n[i,k])
      logit(p[i,k])<-  theta[i,k]
      theta[i,k] <- u[i] + delta[i,k] # Trial specific u + EMaxModel from RE
    }

    for (k in 2:na[i]){
      %s # linkstr
      #delta[i,k] ~ dnorm(d[i,k],prec)
      delta[i,k] <- d[i,k]
    }
  }

  %s # prior for b1
  %s # prior for b2
  # heterogeniety across trials
  tau~dunif(0,5)      # common standard deviation is given a vague prior
  prec<-1/pow(tau,2) # common precision in RE-NMA model

  ",dose.str,prior.b1.str,prior.b2.str)
  }
  jagsmodelDRnetmeta <- (eval(parse(text=paste0('jagsmodelDRnetmeta <- function(){',code.str,'}'))))
return(jagsmodelDRnetmeta)
}


# mym <- makeJAGSmodel(metareg=F,effects='random',consistency=F,class=F,doselink='spline')
#
# smodel <- jags.model(textConnection(mym),data = dosresnetmeta.jagsdata)
#
# #newfunctiontext="mym<-function()"
#
# dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','tau','p.drug','dev','rhat','r','n','totresdev','Z'),model.file = modelDRnetmetaSpline,
#                                           n.chains=3,n.iter = 100,n.burnin = 20,DIC=F,n.thin = 1)
#
# modelDRnetmetaSpline <- (eval(parse(text=paste0('modelDRnetmetaSpline <- function(){',mym,'}'))))
#
# # %s # prior for d2
# # %s # prior for b