getEstimate <- function(model) {
   par_estimates <- parameterEstimates(model, standardized = TRUE)
   par_estimates$newlabel <- ifelse(par_estimates$label == "", 
                                    paste(par_estimates$lhs, par_estimates$op,par_estimates$rhs, sep=""),
                                    par_estimates$label)
   result <- par_estimates$est
   ses <- par_estimates$se
   names(result) <- par_estimates$newlabel
   names(ses) <- par_estimates$newlabel
   return(list(estimates = result[which(!is.na(par_estimates$z))],
               ses = ses[which(!is.na(par_estimates$z))]))
}

MCMed <- function (expression, ..., ACM = NULL, object = NULL, rep = 5000,
          CI = 95, seed = 2022) {
  set.seed(seed)
  input <- c(...)
  paramnames <- all.vars(stats::as.formula(paste("~",
                                                 expression)))
  if (class(object) == "lavaan") {
    input <- getEstimate(object)$estimates[paramnames]
    ses <- getEstimate(object)$ses[paramnames]
    # ACM <- lavaan::vcov(object)[paramnames, paramnames]
  }
  vecs <- list()
  for (i in 1:length(input)) {
    item <- rnorm(rep) * ses[i] + input[i]
    item <- c(input[i], item)
    vecs[[names(input[i])]] <- item
  }

  # dat <- MASS::mvrnorm(n = rep, mu = input, Sigma = ACM)
  # dat <- rbind(input, dat)

  # vecs <- as.list(as.data.frame(dat))
  for (i in 1:length(vecs)) {
    assign(paramnames[i], vecs[[i]])
  }
  indirect <- eval(parse(text = expression))
  low <- (1 - CI/100)/2
  upp <- ((1 - CI/100)/2) + (CI/100)
  LL <- round(quantile(indirect[-1], low), digits = 3)
  UL <- round(quantile(indirect[-1], upp), digits = 3)
  interval <- list(round(indirect[1], 3), rbind(LL, UL))
  dimnames(interval[[2]]) <- list(c("LL", "UL"),
                                  c(" "))
  names(interval) <- c("Point Estimate", paste(CI, "% Confidence Interval",
                                               sep = ""))

  return(interval)
}

# MCMed <- function (expression, ..., ACM = NULL, object = NULL, rep = 20000, 
#           CI = 95) {
#   input <- c(...)
#   paramnames <- all.vars(stats::as.formula(paste("~", 
#                                                  expression)))
#   if (class(object) == "lavaan") {
#     input <- getEstimate(object)[paramnames]
#     ACM <- lavaan::vcov(object)[paramnames, paramnames]
#   }
#   print(input)
#   vecs <- list()
#   dat <- MASS::mvrnorm(n = rep, mu = input, Sigma = ACM)
#   dat <- rbind(input, dat)
#   vecs <- as.list(as.data.frame(dat))
#   for (i in 1:length(vecs)) {
#     assign(paramnames[i], vecs[[i]])
#   }
#   indirect <- eval(parse(text = expression))
#   low <- (1 - CI/100)/2
#   upp <- ((1 - CI/100)/2) + (CI/100)
#   LL <- round(quantile(indirect[-1], low), digits = 4)
#   UL <- round(quantile(indirect[-1], upp), digits = 4)
#   interval <- list(indirect[1], rbind(LL, UL))
#   dimnames(interval[[2]]) <- list(c("LL", "UL"), 
#                                   c(" "))
#   names(interval) <- c("Point Estimate", paste(CI, "% Confidence Interval", 
#                                                sep = ""))
# 
#   return(interval)
# }
