R/summarize_simulation2.R

#' Summarize MCMC Output (for one simulation run)
#'
#'A function that intakes MCMC samples and outputs a tibble that consists of posterior estimates with their fixed values at start in one tibble.
#'Other variable included in tibble is bias calculated for each parameter. Included parameters in the tibble are all global parameters, country-year estimates of sens, spec, vradj and gamma.truemat.vr.
#' @param simsetting Describes sim setting
#' @param runname Describes run setting
#' @param simnum Simulation number i out of Nsim
#' @param globalvars Vector of global variable names
#' @param output.dir Output directory name
#'
#' @return Posterior estimates of global variables and country-year parameters (sens, spec, vradj, gammatruematvr), calculate bias, and set fixed values from simulation code.
#' @export
#'
#' @examples
#' summarize_simulation(simsetting = "test", runname = "test", simnum = 1, globalvars = c("sensworld"), output.dir = output.dir)



summarize_simulation <-
  function(simsetting,
           runname,
           simnum,
           globalvars = globalvars,
           output.dir) {
    tib <-
      readRDS(paste(output.dir, 'mcmctib_', simnum, '.RDS', sep = ""))

    tidydat <- tib %>%
      select(
        contains("sens"),
        contains("spec"),
        contains("vradj.ct"),
        contains("gamma.truematvr.ct")
      )

    if (ncol(tidydat) < 62) { #for testing purposes
      #only and all selected parameter should be in this set
      stop(
        paste0(
          "At least oneparameter of interest is missing. .\n",
          "Check to make sure mcmctib_ is correct"
        )
      )
      call. = FALSE
    } else {
      tidydat <- tidydat %>%
        gather(key = "parameter", value = "value") %>%
        group_by(parameter) %>%
        summarize(
          lower = quantile(value, .025),
          median = quantile(value, .5),
          upper = quantile(value, .975)
        ) %>%
        mutate(sim = simnum, runsetting = runname
        )

    #Pulls in simlist with the set values, but this is a list so much call by name
    truthdat <-
      readRDS(paste(output.dir, "simlist_iter", simnum, ".RDS", sep = ""))
    global_settings <-
      readRDS(paste(output.dir, "global_settings.RDS", sep = ""))

    globalvars <-
      c("sensworld",
        "specworld",
        "rho.alphabeta",
        "sigma.alpha",
        "sigma.beta")


    C <- truthdat$C
    Nyears <- truthdat$Nyears
    tidydat[["truth"]] <- NA

    for (c in 1:C) {
      for (t in 1:Nyears) {
        tidydat[which(tidydat$parameter == paste("gamma.truematvr.ct[", c, ",", t , "]", sep = "")), "truth"] <-  truthdat$GAMMATRUE.ct[c, t]

        tidydat[which(tidydat$parameter == paste("sens.ct[", c, ",", t , "]" , sep = "")), "truth"] <- truthdat$SENS.ct[c, t]

         tidydat[which(tidydat$parameter == paste("spec.ct[", c, ",", t , "]", sep = "")), "truth"] <- truthdat$SPEC[c, t]

        tidydat[which(tidydat$parameter == paste("vradj.ct[", c, ",", t , "]", sep = "")), "truth"] <-  truthdat$VRADJ.ct[c, t]

      }
    }

    for (g in 1:length(globalvars)) {
      tidydat[which(tidydat$parameter == globalvars[g]),  "truth"] <-
        global_settings[[globalvars[g]]]
    }

    tidydat <- tidydat %>%
      mutate(bias = (truth - median),
             coverage = NA)

    for (i in 1:nrow(tidydat)) {
      if ((tidydat$lower[i] < tidydat$truth[i]) == TRUE &&
          (tidydat$upper[i] > tidydat$truth[i]) == TRUE) {
        tidydat$coverage[i] <- 1
      } else {
        tidydat$coverage[i] <- 0
      }
    }


    saveRDS(tidydat, file = paste(output.dir, "sumtib_", simnum, '.RDS', sep = ""))

    }
  }
Enpeterson/outputsim documentation built on May 24, 2019, 9:53 a.m.