outputsim: Output Simulation Results

#' Run Simulation Output
#'
#'Takes in global set values for global parameters, and simulates country-year parameters, which feeds into simulated country-year data.
#'A list is outputted that contains both simulated country-year data to be input for JAGS run, and fixed "truth" country-year parameters values ie SENS.ct.
#'Capitolized parameters in list are indicating fixed values for country-years to not interfere with jags runs.
#'
#' @param runname Description of run settings
#' @param i Simulation number out of Nsim
#' @param global_settings a vector of fixed values for global parameters, ie sensworld, specworld
#' @param simsetting Description of simsetting
#'
#' @return returns a list which consists of simulated data and jags inputs.
#' @export
#'
#' @examples
#' run_simulation( runname = "test", i=1, global_settings = global_settings, simsetting = "test)




run_simulation <- function(runname , i,
                          global_settings,
                          simsetting){


simlist <- list()
      sensworld <- global_settings[["sensworld"]]
      specworld <- global_settings[["specworld"]]
      rsigma.alpha <- global_settings[["rsigma.alpha"]]
      rsigma.beta <- global_settings[["rsigma.beta"]]
      sigma.alpha <- global_settings[["sigma.alpha"]]
      sigma.beta <- global_settings[["sigma.beta"]]
      rho.alphabeta <- global_settings[["rho.alphabeta"]]
      rrho.alphabeta <- global_settings[["rrho.alphabeta"]]
      tref <- global_settings[["refyear"]]
      C = global_settings[["C"]]
      Nyears = global_settings[["Nyears"]]
      vrenv = global_settings[["vrenv"]]
      K=2
      B=4

      covart <- rho.alphabeta * sigma.beta * sigma.alpha
      condsd_spec <- sqrt(sigma.beta^2 - (covart^2/sigma.alpha^2))
      gammatrue.ct <- cond_mean_lambda.ct <-  array(NA, c(C, Nyears))
      lambda.ctk <-  mean_lambda.ctk <-    array(NA, c(C, Nyears, K))


      #Simulate parameter values

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

          set.seed(1234)


          #Set gammatrue to G for all years
         gammatrue.ct[c,t] <- 0.01 #Set gammatrue to G for all years
        }
         #Set sens and spec for ref year = global sens and spec with no var because only simming one country.
         lambda.ctk[c,tref,1] <- rtruncnorm(1, a=0, b=0.7, mean = sensworld,  sd = rsigma.alpha)
         lambda.ctk[c,tref,2] <- rtruncnorm(1, a=0, b=1, mean = specworld ,  sd = rsigma.beta)

         #Create draws from conditional RW for sens and spec
         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]<- rtruncnorm(1, a = 0.1, b = 0.7, mean = mean_lambda.ctk[c,tindex,1],  sd = sigma.alpha)
           cond_mean_lambda.ct[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] <-  rtruncnorm(1, a = 0.8, b = 1, cond_mean_lambda.ct[c,tindex],  condsd_spec)
         }
         for(tindex in (tref-1):1){
           mean_lambda.ctk[c,tindex,1:K] <- lambda.ctk[c,tindex + 1,1:K]
           lambda.ctk[c,tindex,1]<- rtruncnorm(1, a = 0.1, b = 0.7, mean = mean_lambda.ctk[c,tindex,1],  sd = sigma.alpha)
           cond_mean_lambda.ct[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] <-  rtruncnorm(1, a = 0.8, b = 1, cond_mean_lambda.ct[c,tindex],  condsd_spec)
         }

      }



    #Generate multinomial data and cumulative counts
      matvr.ct <- vr.ct <- TP.ct <- FP.ct <- FN.ct <- TN.ct <- rhovr.ct <- truematvr.ct <- array(NA, c(C, Nyears))
      rho.ctb <- array(NA, c(C, Nyears, B))
      SENS.ct <- SPEC.ct <- VRADJ.ct <- array(NA, c(C, Nyears))

      S <- 4
      tp <- 1
      fp <- 2
      tn <- 3
      fn <- 4


      for(c in 1:C){
        for(t in 1:Nyears){
          rhovr.ct[c,t] <- 1
          rho.ctb[c,t,tp] =  lambda.ctk[c,t,1] * gammatrue.ct[c,t]* rhovr.ct[c,t]
          rho.ctb[c,t,fn] = gammatrue.ct[c,t]* rhovr.ct[c,t] - rho.ctb[c,t,tp]
          rho.ctb[c,t,tn]  = lambda.ctk[c,t,2]* (1-gammatrue.ct[c,t])* rhovr.ct[c,t]
          rho.ctb[c,t,fp] = (1-gammatrue.ct[c,t])* rhovr.ct[c,t] - rho.ctb[c,t,tn]
          SENS.ct[c,t] <- lambda.ctk[c,t,1]
          SPEC.ct[c,t] <- lambda.ctk[c,t,2]
          VRADJ.ct[c,t] <- sum(rho.ctb[c,t,c(tp, fn)])/ sum(rho.ctb[c,t,c(tp, fp)])
          }
      }
      multinominvr.cts <- array(NA, c(C, Nyears, S))
      set.seed(i)
      for(c in 1:C){
        for(t in 1:Nyears){
          vr.ct[c,t] <- rbinom(1, vrenv, prob = runif(1, 0.2,1))
          multinominvr.cts[c,t,1:S] <- rmultinom(1, vr.ct[c,t] , prob = c(rho.ctb[c,t,1:S]) )
          matvr.ct[c,t] <- sum(multinominvr.cts[c,t, c(tp, fp)])
          truematvr.ct[c,t] <- sum(multinominvr.cts[c,t, c(tp, fn)])
          TP.ct[c,t] <- multinominvr.cts[c,t,tp]
          TN.ct[c,t] <- multinominvr.cts[c,t,tn]
          FP.ct[c,t] <- multinominvr.cts[c,t, fp]
          FN.ct[c,t] <- multinominvr.cts[c,t, fn]
        }
      }



simlist[["TP.ct"]] <- TP.ct
simlist[["FP.ct"]] <- FP.ct
simlist[["FN.ct"]] <- FN.ct
simlist[["TN.ct"]] <- TN.ct
simlist[["S"]] <- S
simlist[["SENS.ct"]] <- SENS.ct
simlist[["SPEC.ct"]] <- SPEC.ct
simlist[["VRADJ.ct"]] <- VRADJ.ct
simlist[["GAMMATRUE.ct"]] <- gammatrue.ct

if(runname %in% c(simsetting)){
simlist[["multinominvr.cts"]] <- multinominvr.cts
}
simlist[["C"]] <- C
simlist[["Nyears"]] <- Nyears
simlist[["tref"]] <- tref
simlist[["K"]] <- K
simlist[["matvr.ct"]] <- matvr.ct
simlist[["truematvr.ct"]] <- truematvr.ct

simlist[["vr.ct"]] <- vr.ct
simlist[["B"]] <- B
simlist[["tp"]] <- tp
simlist[["fp"]] <- fp
simlist[["tn"]] <- tn
simlist[["fn"]] <- fn
simlist[["indices.matvr"]] <- c(tp, fp)
simlist[["indices.truematvr"]] <- c(tp, fn)

saveRDS(simlist, paste(output.dir, "simlist_iter",i,  ".RDS", sep=""))
return(NULL)
}#end function