R/syn_ordinal.R

In RyanHornby/LipschitzWeights: What the Package Does (One Line, Title Case)

#'
#' @import rstanarm
#' @import gridExtra
#' @import tidyverse
#' @import reshape2
#' @export

syn_ordinal = function(origonal_data, model_string = "outcome ~ 1", chains = 1, iterations = 1000,
                        scale_factor = 1.0, shift_factor = 0.0, c = 0.95, thresh = 1.00,
                        risks_method = "linear", m = 20, thin = 5, plots = FALSE) {
  
  require(rstanarm)
  
  fit_notw <- stan_glm(
    model_string,
    data = origonal_data,
    method = "probit",
    prior = R2(0.2, "mean"), 
    refresh = 0,
    chains = chains, iter = iterations
  )
  
  syn1 = syn_part_1(fit_notw, scale_factor, shift_factor)
  
  weights_LW <- syn1[[2]]
  
  assign("weights_LW", weights_LW, sys.nframe())
  
  fit_LW <- stan_glm(
    model_string,
    data = origonal_data,
    method = "probit",
    prior = R2(0.2, "mean"), 
    weights = weights_LW,
    refresh = 0,
    chains = chains, iter = iterations
  )
  
  rm(list = c("weights_LW"), envir = globalenv())
  
  syn2 = syn_part_2(fit_LW, scale_factor, shift_factor, weights_LW, c)
  
  weights_LW_final = syn2[[2]]
  
  assign("weights_LW_final", weights_LW_final, globalenv())
  
  fit_LW_final <- stan_glm(
    model_string,
    data = origonal_data,
    method = "probit",
    prior = R2(0.2, "mean"), 
    weights = weights_LW_final,
    refresh = 0,
    chains = chains, iter = iterations
  )
  
  rm(list = c("weights_LW_final"), envir = globalenv())
  
  syn3 = syn_part_3(fit_LW_final, scale_factor, shift_factor)
  
  #### synthesis ####
  N = length(origonal_data[,1])
  draws <- as.data.frame(fit_LW_final)
  start <- 500 / 2
  syndata <- vector("list", m)
  for (i in 1:m){
    draws_exp <- exp(draws[start + thin * (i - 1), ])
    syndata[[i]] <- rpois(N, draws_exp)
  }
  
  
  if (plots) {
    syn_plots(syn1[[1]], syn2[[1]], syn3[[1]], syn1[[2]], syn2[[2]])
  }
  
  return(syndata)
}