quick_sim_slopes <- function(model, pred, modx, modx.values) {
  library(lme4)
  sumoutmodel <- summary(model)
  psi <- coef(sumoutmodel)[,1]
  vpsi <- vcov(sumoutmodel)
  psiname <- names(psi)
  
  cint <- matrix(0, length(psi), length(modx.values))
  cslope <- matrix(0, length(psi), length(modx.values))
  
  cintpos1 <- which(psiname == "(Intercept)")
  cintposvec <- which(psiname == modx)
  cslopepos1 <- which(psiname == pred)
  cslopeposvec <- c(which(psiname == paste0(pred, ":", modx)), which(psiname == paste0(modx, ":", pred)))
  
  cint[cintpos1, ] <- 1
  cint[cintposvec,] <- modx.values
  cslope[cslopepos1, ] <- 1
  cslope[cslopeposvec,] <- modx.values
  
  cintest <- t(cint) %*% psi
  cintse <- sqrt(diag(t(cint) %*% vpsi %*% cint))
  cintzval <- cintest/cintse
  cintpval <- round(pnorm(-abs(cintzval)) * 2, 6)
  cintresult <- data.frame(modx.values, est=cintest, se=cintse, z=cintzval, p=cintpval)
  
  cslopeest <- t(cslope) %*% psi
  cslopese <- sqrt(diag(t(cslope) %*% vpsi %*% cslope))
  cslopezval <- cslopeest/cslopese
  cslopepval <- round(pnorm(-abs(cslopezval)) * 2, 6)
  csloperesult <- data.frame(modx.values, est=cslopeest, se=cslopese, z=cslopezval, p=cslopepval)
  
  result <- list("Conditional Intercept" = cintresult, "Conditional Slope" = csloperesult)
  return(result)
}