outputsim: Output Simulation Results

#' Write Bayesian Model
#'A fuction that writes a Bayesian model for sensitivity and specificity based on given settings.
#'
#' @param output.dir The output directory where the model.txt file will be saved.
#' @param modelonegammatrue A TRUE or FALSE indicator, if TRUE the Bayesian model includes one parameter for gammatrue.
#' @param modelgammatruect A TRUE or FALSE indicator, if TRUE the Bayesian model includes a unique gammatrue for each country-year.
#' @param simsetting Saves write model file under the simsetting name.
#'
#' @return a model.txt file that specifies process and data models, and is called within funciton run_mcmc. The model.txt file is saved within the output folder.
#' @export
#'
#' @examples write_model_to_run(output.dir, modelonegammatrue = FALSE, modelgammatruect=FALSE, simsetting = "test")





write_model_to_run <- function(output.dir, modelonegammatrue = FALSE, modelgammatruect=FALSE, simsetting ){


JAGSmodel.name <- "run_data_model.txt"

  cat("

      model{",
      sep = "",
      append = FALSE,
      file = file.path(output.dir, JAGSmodel.name),
      fill = TRUE
  )

  if (modelonegammatrue) {
    cat("
    gamma_true ~ dunif(0, 1)

    for(c in 1:C){
    for(t in 1:Nyears){
    gamma.truematvr.ct[c,t] <- gamma_true
    }
    }
", sep = "", append = TRUE, file = file.path(output.dir, JAGSmodel.name), fill = TRUE)
  }

  if(modelgammatruect) {
    cat("
        for(c in 1:C){
        for(t i 1:Nyears){
        gamma.truematvr.ct[c,t] ~ dunif(0, 1)
  }
}
        ", sep = "", append = TRUE, file = file.path(output.dir, JAGSmodel.name), fill = TRUE)
  }

cat("
for(c in 1:C){
 for(t in 1:Nyears){
    sens.ct[c,t] <- lambda.ctk[c,t,1]
    spec.ct[c,t] <- lambda.ctk[c,t,2]
    }


# #Using MVN conditioning spec on sens
# mean_lambda.ctk[c,tref,1] <- sensworld
# mean_lambda.ctk[c,tref,2] <- specworld

lambda.ctk[c,tref,1] <- sensworld
lambda.ctk[c,tref,2] <- specworld

for(tindex in (tref+1):Nyears){
mean_lambda.ctk[c,tindex,1:K] <- lambda.ctk[c,tindex-1,1:K]
lambda.ctk[c,tindex,1] ~ dnorm(mean_lambda.ctk[c,tindex,1],  prec_eta1)T(0.1,1 )
cond_mean_lambda.ctk[c,tindex] <- mean_lambda.ctk[c,tindex,2] + sigma.beta/sigma.alpha*rho.alphabeta*(lambda.ctk[c,tindex,1] - mean_lambda.ctk[c,tindex,1])
lambda.ctk[c,tindex,2] ~ dnorm(cond_mean_lambda.ctk[c,tindex],  condprec_eta2)T(0.90,1)
}

for(tindex in 2:tref){
mean_lambda.ctk[c,tindex-1,1:K] <- lambda.ctk[c,tindex,1:K]
lambda.ctk[c,tindex-1,1] ~ dnorm(mean_lambda.ctk[c,tindex-1,1],  prec_eta1)T(0.1,1)
cond_mean_lambda.ctk[c,tindex-1] <- mean_lambda.ctk[c,tindex-1,2] + sigma.beta/sigma.alpha*rho.alphabeta*(lambda.ctk[c,tindex-1,1] - mean_lambda.ctk[c,tindex-1,1])
lambda.ctk[c,tindex-1,2] ~ dnorm(cond_mean_lambda.ctk[c,tindex-1],  condprec_eta2)T(0.95,1)
}


} #end C loop

#prec_eta1tref <- 1/rsigma.alpha^2
prec_eta1 <- 1/sigma.alpha^2
#prec_eta2tref <- 1/rsigma.beta^2
prec_eta2 <- 1/sigma.beta^2



#condprec_eta2tref <- 1/((1-rrho.alphabeta^2)*rsigma.beta^2)
#covartref <- rrho.alphabeta * rsigma.beta * rsigma.alpha
#condprec_eta2tref <- 1/(rsigma.beta^2 - (covartref^2/rsigma.alpha^2))

#condprec_eta2 <- 1/((1-rho.alphabeta^2)*sigma.beta^2)
covart <- rho.alphabeta * sigma.beta * sigma.alpha
condprec_eta2 <- 1/(sigma.beta^2 - (covart^2/sigma.alpha^2))

#Sigma_ref[1,1] <- rsigma.alpha^2
#Sigma_ref[2,2] <- rsigma.beta^2
#Sigma_ref[1,2] <- rho.alphabeta * rsigma.alpha * rsigma.beta
#Sigma_ref[2,1] <- Sigma_ref[1,2]


Sigma_rw[1,1] <- sigma.alpha^2
Sigma_rw[2,2] <- sigma.beta^2
Sigma_rw[1,2] <- rho.alphabeta * sigma.alpha * sigma.beta
Sigma_rw[2,1] <- Sigma_rw[1,2]

sigma.alpha ~ dunif(0,5)
#rsigma.alpha ~ dunif(0,5)
rsigma.alpha <- 0

sigma.beta ~ dunif(0,5)
#rsigma.beta ~ dunif(0,5)
rsigma.beta <- 0

rho.alphabeta ~ dunif(-1,1)
rrho.alphabeta <- 0

sensworld ~ dunif(0,1)
specworld ~ dunif(0,1)

", sep = "", append = TRUE, file = file.path(output.dir, JAGSmodel.name), fill = TRUE)





#--------------------- Reparametrize
#----------------------------------------------------------------------------------------------------------------------


cat("

    for (c in 1:C){
    for (t in 1:Nyears){

    # also calculate vr adjustment
    # based on estimated gammamatinvr (based on observed may be an issue so use function in R)
    vradj.ct[c,t] <- (rho.ctb[c,t,fn] + rho.ctb[c,t,tp])/gamma.matvr.ct[c,t]


    gamma.matvr.ct[c,t] <- sens.ct[c,t]*gamma.truematvr.ct[c,t] + (1-spec.ct[c,t])*(1-gamma.truematvr.ct[c,t])

    rhovr.ct[c,t] <- 1
    rho.ctb[c,t,tp] =  sens.ct[c,t] * gamma.truematvr.ct[c,t]* rhovr.ct[c,t]
    rho.ctb[c,t,fn] = gamma.truematvr.ct[c,t]* rhovr.ct[c,t] - rho.ctb[c,t,tp]
    rho.ctb[c,t,tn]  = spec.ct[c,t]* (1-gamma.truematvr.ct[c,t])* rhovr.ct[c,t]
    rho.ctb[c,t,fp] = (1-gamma.truematvr.ct[c,t])* rhovr.ct[c,t] - rho.ctb[c,t,tn]

    }
}

    ",
    sep = "",
    append = TRUE,
    file = file.path(output.dir, JAGSmodel.name),
    fill = TRUE
)



#--------------------- Var-Cov matrix
#----------------------------------------------------------------------------------------------------------------------


cat("


   # define the varcov matrix for the multnomial for all c,t

    for(c in 1:C){
    for(t in 1:Nyears){
    for (b in 1:B){

    Sigma.single.ctbb[c,t,b,b] <- rho.ctb[c,t,b]*(1-rho.ctb[c,t,b])/vr.ct[c,t]

    }
    for (b1 in 1:(B-1)){
    for (b2 in (b1+1):B){

    Sigma.single.ctbb[c,t,b1,b2] <- -(rho.ctb[c,t,b1]*rho.ctb[c,t,b2])/vr.ct[c,t]
    Sigma.single.ctbb[c,t,b2,b1] <- Sigma.single.ctbb[c,t,b1,b2]

    }} #end b loop
    }#end t loop
    } #end c loop
 ",
    sep = "",
    append = TRUE,
    file = file.path(output.dir, JAGSmodel.name),
    fill = TRUE
)



#--------------------- data models
#----------------------------------------------------------------------------------------------------------------------

if(simsetting %in% c("alldata")) {
cat("

for(c in 1:C){
for(t in 1:Nyears){

    multinominvr.cts[c,t,1:S] ~ dmulti(mumultiinvr.cts[c,t,1:S], vr.ct[c,t])
    mumultiinvr.cts[c,t,tp] <- rho.ctb[c, t, tp]
    mumultiinvr.cts[c,t,fp] <- rho.ctb[c, t, fp]
    mumultiinvr.cts[c,t,tn] <- rho.ctb[c, t, tn]
    mumultiinvr.cts[c,t,fn] <- rho.ctb[c, t, fn]
    }#end t loop
}#end C loop




    ",
    sep = "",
    append = TRUE,
    file = file.path(output.dir, JAGSmodel.name),
    fill = TRUE
)
}
if(simsetting %in% c("nobreakdown")) {
  cat("

      for(c in 1:C){
      for(t in 1:Nyears){
      matvr.ct[c,t] ~ dbinom(rhomatvr.ct[c,t], vr.ct[c,t])
      rhomatvr.ct[c,t] <- sum(rho.ctb[c,t,indices.matvr])

      truematvr.ct[c,t] ~ dbinom(rhotruematvr.ct[c,t], vr.ct[c,t])
      rhotruematvr.ct[c,t] <- sum(rho.ctb[c,t,indices.truematvr])

      }#end t loop
      }#end C loop




      ",
      sep = "",
      append = TRUE,
      file = file.path(output.dir, JAGSmodel.name),
      fill = TRUE
  )
}




  cat("
      #   ----------------------------------
}  #end data




      ",
      sep = "",
      append = TRUE,
      file = file.path(output.dir, JAGSmodel.name),
      fill = TRUE
  )
  return(NULL)
}