inst/shiny/skprGUI/server.R

In skpr: Design of Experiments Suite: Generate and Evaluate Optimal Designs

library(skpr)
library(rintrojs)
library(gt)

function(input, output, session) {

  inc_progress_session = function(amount = 0.1, message = NULL, detail = NULL) incProgress(amount, message, detail, session)

  inputlist_htc = reactive({
    input$submitbutton
    inputlist1 = list()
    for (i in 1:6) {
      if ( (i == 1 && input$numberfactors > 0) && isolate(input$blockdepth1) == "htc") {
        if (input$factortype1 == "numeric") {
          inputlist1[[input$factorname1]] = seq(input$numericlow1, input$numerichigh1, length.out = input$numericlength1)
        }
        if (input$factortype1 == "discnum") {
          inputlist1[[input$factorname1]] = as.numeric(strsplit(gsub(" ", "", input$disclevels1, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype1 == "cat") {
          inputlist1[[input$factorname1]] = strsplit(gsub(" ", "", input$levels1, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 2 && input$numberfactors > 1) && isolate(input$blockdepth2) == "htc") {
        if (input$factortype2 == "numeric") {
          inputlist1[[input$factorname2]] = seq(input$numericlow2, input$numerichigh2, length.out = input$numericlength2)
        }
        if (input$factortype2 == "discnum") {
          inputlist1[[input$factorname2]] = as.numeric(strsplit(gsub(" ", "", input$disclevels2, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype2 == "cat") {
          inputlist1[[input$factorname2]] = strsplit(gsub(" ", "", input$levels2, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 3 && input$numberfactors > 2) && isolate(input$blockdepth3) == "htc") {
        if (input$factortype3 == "numeric") {
          inputlist1[[input$factorname3]] = seq(input$numericlow3, input$numerichigh3, length.out = input$numericlength3)
        }
        if (input$factortype3 == "discnum") {
          inputlist1[[input$factorname3]] = as.numeric(strsplit(gsub(" ", "", input$disclevels3, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype3 == "cat") {
          inputlist1[[input$factorname3]] = strsplit(gsub(" ", "", input$levels3, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 4 && input$numberfactors > 3) && isolate(input$blockdepth4) == "htc") {
        if (input$factortype4 == "numeric") {
          inputlist1[[input$factorname4]] = seq(input$numericlow4, input$numerichigh4, length.out = input$numericlength4)
        }
        if (input$factortype4 == "discnum") {
          inputlist1[[input$factorname4]] = as.numeric(strsplit(gsub(" ", "", input$disclevels4, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype4 == "cat") {
          inputlist1[[input$factorname4]] = strsplit(gsub(" ", "", input$levels4, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 5 && input$numberfactors > 4) && isolate(input$blockdepth5) == "htc") {
        if (input$factortype5 == "numeric") {
          inputlist1[[input$factorname5]] = seq(input$numericlow5, input$numerichigh5, length.out = input$numericlength5)
        }
        if (input$factortype5 == "discnum") {
          inputlist1[[input$factorname5]] = as.numeric(strsplit(gsub(" ", "", input$disclevels5, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype5 == "cat") {
          inputlist1[[input$factorname5]] = strsplit(gsub(" ", "", input$levels5, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 6 && input$numberfactors > 5) && isolate(input$blockdepth6) == "htc") {
        if (input$factortype6 == "numeric") {
          inputlist1[[input$factorname6]] = seq(input$numericlow6, input$numerichigh6, length.out = input$numericlength6)
        }
        if (input$factortype6 == "discnum") {
          inputlist1[[input$factorname6]] = as.numeric(strsplit(gsub(" ", "", input$disclevels6, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype6 == "cat") {
          inputlist1[[input$factorname6]] = strsplit(gsub(" ", "", input$levels6, fixed = TRUE), split = ",")[[1]]
        }
      }
    }
    inputlist1
  })

  inputlist_htctext = reactive({
    inputlist1 = list()
    for (i in 1:6) {
      if ( (i == 1 && input$numberfactors > 0) && input$blockdepth1 == "htc") {
        if (input$factortype1 == "numeric") {
          inputlist1[[input$factorname1]] = seq(input$numericlow1, input$numerichigh1, length.out = input$numericlength1)
        }
        if (input$factortype1 == "discnum") {
          inputlist1[[input$factorname1]] = as.numeric(strsplit(gsub(" ", "", input$disclevels1, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype1 == "cat") {
          inputlist1[[input$factorname1]] = strsplit(gsub(" ", "", input$levels1, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 2 && input$numberfactors > 1) && input$blockdepth2 == "htc") {
        if (input$factortype2 == "numeric") {
          inputlist1[[input$factorname2]] = seq(input$numericlow2, input$numerichigh2, length.out = input$numericlength2)
        }
        if (input$factortype2 == "discnum") {
          inputlist1[[input$factorname2]] = as.numeric(strsplit(gsub(" ", "", input$disclevels2, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype2 == "cat") {
          inputlist1[[input$factorname2]] = strsplit(gsub(" ", "", input$levels2, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 3 && input$numberfactors > 2) && input$blockdepth3 == "htc") {
        if (input$factortype3 == "numeric") {
          inputlist1[[input$factorname3]] = seq(input$numericlow3, input$numerichigh3, length.out = input$numericlength3)
        }
        if (input$factortype3 == "discnum") {
          inputlist1[[input$factorname3]] = as.numeric(strsplit(gsub(" ", "", input$disclevels3, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype3 == "cat") {
          inputlist1[[input$factorname3]] = strsplit(gsub(" ", "", input$levels3, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 4 && input$numberfactors > 3) && input$blockdepth4 == "htc") {
        if (input$factortype4 == "numeric") {
          inputlist1[[input$factorname4]] = seq(input$numericlow4, input$numerichigh4, length.out = input$numericlength4)
        }
        if (input$factortype4 == "discnum") {
          inputlist1[[input$factorname4]] = as.numeric(strsplit(gsub(" ", "", input$disclevels4, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype4 == "cat") {
          inputlist1[[input$factorname4]] = strsplit(gsub(" ", "", input$levels4, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 5 && input$numberfactors > 4) && input$blockdepth5 == "htc") {
        if (input$factortype5 == "numeric") {
          inputlist1[[input$factorname5]] = seq(input$numericlow5, input$numerichigh5, length.out = input$numericlength5)
        }
        if (input$factortype5 == "discnum") {
          inputlist1[[input$factorname5]] = as.numeric(strsplit(gsub(" ", "", input$disclevels5, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype5 == "cat") {
          inputlist1[[input$factorname5]] = strsplit(gsub(" ", "", input$levels5, fixed = TRUE), split = ",")[[1]]
        }
      }
      if ( (i == 6 && input$numberfactors > 5) && input$blockdepth6 == "htc") {
        if (input$factortype6 == "numeric") {
          inputlist1[[input$factorname6]] = seq(input$numericlow6, input$numerichigh6, length.out = input$numericlength6)
        }
        if (input$factortype6 == "discnum") {
          inputlist1[[input$factorname6]] = as.numeric(strsplit(gsub(" ", "", input$disclevels6, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype6 == "cat") {
          inputlist1[[input$factorname6]] = strsplit(gsub(" ", "", input$levels6, fixed = TRUE), split = ",")[[1]]
        }
      }
    }
    inputlist1
  })


  inputlist = reactive({
    inputlist1 = list()
    for (i in 1:6) {
      if (i == 1 && input$numberfactors > 0 ) {
        if (input$blockdepth1 == "etc") {
          if (input$factortype1 == "numeric") {
            inputlist1[[input$factorname1]] = seq(input$numericlow1, input$numerichigh1, length.out = input$numericlength1)
          }
          if (input$factortype1 == "discnum") {
            inputlist1[[input$factorname1]] = as.numeric(strsplit(gsub(" ", "", input$disclevels1, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype1 == "cat") {
            inputlist1[[input$factorname1]] = strsplit(gsub(" ", "", input$levels1, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
      if (i == 2 && input$numberfactors > 1) {
        if (input$blockdepth2 == "etc") {
          if (input$factortype2 == "numeric") {
            inputlist1[[input$factorname2]] = seq(input$numericlow2, input$numerichigh2, length.out = input$numericlength2)
          }
          if (input$factortype2 == "discnum") {
            inputlist1[[input$factorname2]] = as.numeric(strsplit(gsub(" ", "", input$disclevels2, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype2 == "cat") {
            inputlist1[[input$factorname2]] = strsplit(gsub(" ", "", input$levels2, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
      if (i == 3 && input$numberfactors > 2) {
        if (input$blockdepth3 == "etc") {
          if (input$factortype3 == "numeric") {
            inputlist1[[input$factorname3]] = seq(input$numericlow3, input$numerichigh3, length.out = input$numericlength3)
          }
          if (input$factortype3 == "discnum") {
            inputlist1[[input$factorname3]] = as.numeric(strsplit(gsub(" ", "", input$disclevels3, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype3 == "cat") {
            inputlist1[[input$factorname3]] = strsplit(gsub(" ", "", input$levels3, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
      if (i == 4 && input$numberfactors > 3) {
        if (input$blockdepth4 == "etc") {
          if (input$factortype4 == "numeric") {
            inputlist1[[input$factorname4]] = seq(input$numericlow4, input$numerichigh4, length.out = input$numericlength4)
          }
          if (input$factortype4 == "discnum") {
            inputlist1[[input$factorname4]] = as.numeric(strsplit(gsub(" ", "", input$disclevels4, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype4 == "cat") {
            inputlist1[[input$factorname4]] = strsplit(gsub(" ", "", input$levels4, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
      if (i == 5 && input$numberfactors > 4) {
        if (input$blockdepth5 == "etc") {
          if (input$factortype5 == "numeric") {
            inputlist1[[input$factorname5]] = seq(input$numericlow5, input$numerichigh5, length.out = input$numericlength5)
          }
          if (input$factortype5 == "discnum") {
            inputlist1[[input$factorname5]] = as.numeric(strsplit(gsub(" ", "", input$disclevels5, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype5 == "cat") {
            inputlist1[[input$factorname5]] = strsplit(gsub(" ", "", input$levels5, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
      if (i == 6 && input$numberfactors > 5) {
        if (input$blockdepth6 == "etc") {
          if (input$factortype6 == "numeric") {
            inputlist1[[input$factorname6]] = seq(input$numericlow6, input$numerichigh6, length.out = input$numericlength6)
          }
          if (input$factortype6 == "discnum") {
            inputlist1[[input$factorname6]] = as.numeric(strsplit(gsub(" ", "", input$disclevels6, fixed = TRUE), split = ",")[[1]])
          }
          if (input$factortype6 == "cat") {
            inputlist1[[input$factorname6]] = strsplit(gsub(" ", "", input$levels6, fixed = TRUE), split = ",")[[1]]
          }
        }
      }
    }
    inputlist1
  })

  candidatesetall = reactive({
    candidateset1 = list()
    for (i in 1:6) {
      if (i == 1 && input$numberfactors > 0 ) {
        if (input$factortype1 == "numeric") {
          candidateset1[[input$factorname1]] = seq(input$numericlow1, input$numerichigh1, length.out = input$numericlength1)
        }
        if (input$factortype1 == "discnum") {
          candidateset1[[input$factorname1]] = as.numeric(strsplit(gsub(" ", "", input$disclevels1, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype1 == "cat") {
          candidateset1[[input$factorname1]] = strsplit(gsub(" ", "", input$levels1, fixed = TRUE), split = ",")[[1]]
        }
      }
      if (i == 2 && input$numberfactors > 1) {
        if (input$factortype2 == "numeric") {
          candidateset1[[input$factorname2]] = seq(input$numericlow2, input$numerichigh2, length.out = input$numericlength2)
        }
        if (input$factortype2 == "discnum") {
          candidateset1[[input$factorname2]] = as.numeric(strsplit(gsub(" ", "", input$disclevels2, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype2 == "cat") {
          candidateset1[[input$factorname2]] = strsplit(gsub(" ", "", input$levels2, fixed = TRUE), split = ",")[[1]]
        }
      }
      if (i == 3 && input$numberfactors > 2) {
        if (input$factortype3 == "numeric") {
          candidateset1[[input$factorname3]] = seq(input$numericlow3, input$numerichigh3, length.out = input$numericlength3)
        }
        if (input$factortype3 == "discnum") {
          candidateset1[[input$factorname3]] = as.numeric(strsplit(gsub(" ", "", input$disclevels3, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype3 == "cat") {
          candidateset1[[input$factorname3]] = strsplit(gsub(" ", "", input$levels3, fixed = TRUE), split = ",")[[1]]
        }
      }
      if (i == 4 && input$numberfactors > 3) {
        if (input$factortype4 == "numeric") {
          candidateset1[[input$factorname4]] = seq(input$numericlow4, input$numerichigh4, length.out = input$numericlength4)
        }
        if (input$factortype4 == "discnum") {
          candidateset1[[input$factorname4]] = as.numeric(strsplit(gsub(" ", "", input$disclevels4, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype4 == "cat") {
          candidateset1[[input$factorname4]] = strsplit(gsub(" ", "", input$levels4, fixed = TRUE), split = ",")[[1]]
        }
      }
      if (i == 5 && input$numberfactors > 4) {
        if (input$factortype5 == "numeric") {
          candidateset1[[input$factorname5]] = seq(input$numericlow5, input$numerichigh5, length.out = input$numericlength5)
        }
        if (input$factortype5 == "discnum") {
          candidateset1[[input$factorname5]] = as.numeric(strsplit(gsub(" ", "", input$disclevels5, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype5 == "cat") {
          candidateset1[[input$factorname5]] = strsplit(gsub(" ", "", input$levels5, fixed = TRUE), split = ",")[[1]]
        }
      }
      if (i == 6 && input$numberfactors > 5) {
        if (input$factortype6 == "numeric") {
          candidateset1[[input$factorname6]] = seq(input$numericlow6, input$numerichigh6, length.out = input$numericlength6)
        }
        if (input$factortype6 == "discnum") {
          candidateset1[[input$factorname6]] = as.numeric(strsplit(gsub(" ", "", input$disclevels6, fixed = TRUE), split = ",")[[1]])
        }
        if (input$factortype6 == "cat") {
          candidateset1[[input$factorname6]] = strsplit(gsub(" ", "", input$levels6, fixed = TRUE), split = ",")[[1]]
        }
      }
    }
    candidateset1
  })



  inputstring = reactive({
    updatevector = c(input$blockdepth1, input$blockdepth2, input$blockdepth3, input$blockdepth4, input$blockdepth5, input$blockdepth6)
    commacount = input$numberfactors - 1
    finalstring = c()
    for (i in 1:6) {
      if (i == 1 && input$numberfactors > 0) {
        finalstring = c(finalstring, input$factorname1, " = ")
        if (input$factortype1 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow1, ",", input$numerichigh1, ", length.out = ", input$numericlength1, ")")
        }
        if (input$factortype1 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels1, fixed = TRUE), ")")
        }
        if (input$factortype1 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels1, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
        if (commacount > 0) {
          commacount = commacount - 1
          finalstring = c(finalstring, paste0(c(", \n", rep(" ", 27)), collapse = ""))
        }
      }
      if ( (i == 2 && input$numberfactors > 1)) {
        finalstring = c(finalstring, input$factorname2, " = ")
        if (input$factortype2 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow2, ",", input$numerichigh2, ", length.out = ", input$numericlength2, ")")
        }
        if (input$factortype2 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels2, fixed = TRUE), ")")
        }
        if (input$factortype2 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels2, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
        if (commacount > 0) {
          commacount = commacount - 1
          finalstring = c(finalstring, paste0(c(", \n", rep(" ", 27)), collapse = ""))
        }
      }
      if (i == 3 && input$numberfactors > 2) {
        finalstring = c(finalstring, input$factorname3, " = ")
        if (input$factortype3 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow3, ",", input$numerichigh3, ", length.out = ", input$numericlength3, ")")
        }
        if (input$factortype3 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels3, fixed = TRUE), ")")
        }
        if (input$factortype3 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels3, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
        if (commacount > 0) {
          commacount = commacount - 1
          finalstring = c(finalstring, paste0(c(", \n", rep(" ", 27)), collapse = ""))
        }
      }
      if (i == 4 && input$numberfactors > 3) {
        finalstring = c(finalstring, input$factorname4, " = ")
        if (input$factortype4 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow4, ",", input$numerichigh4, ", length.out = ", input$numericlength4, ")")
        }
        if (input$factortype4 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels4, fixed = TRUE), ")")
        }
        if (input$factortype4 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels4, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
        if (commacount > 0) {
          commacount = commacount - 1
          finalstring = c(finalstring, paste0(c(", \n", rep(" ", 27)), collapse = ""))
        }
      }
      if (i == 5 && input$numberfactors > 4) {
        finalstring = c(finalstring, input$factorname5, " = ")
        if (input$factortype5 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow5, ",", input$numerichigh5, ", length.out = ", input$numericlength5, ")")
        }
        if (input$factortype5 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels5, fixed = TRUE), ")")
        }
        if (input$factortype5 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels5, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
        if (commacount > 0) {
          commacount = commacount - 1
          finalstring = c(finalstring, paste0(c(", \n", rep(" ", 27)), collapse = ""))
        }
      }
      if (i == 6 && input$numberfactors > 5) {
        finalstring = c(finalstring, input$factorname6, " = ")
        if (input$factortype6 == "numeric") {
          finalstring = c(finalstring, "seq(", input$numericlow6, ",", input$numerichigh6, ", length.out = ", input$numericlength6, ")")
        }
        if (input$factortype6 == "discnum") {
          finalstring = c(finalstring, "c(", gsub(" ", "", input$disclevels6, fixed = TRUE), ")")
        }
        if (input$factortype6 == "cat") {
          levelstring = paste0(c("\""), strsplit(gsub(" ", "", input$levels6, fixed = TRUE), split = ",")[[1]], c("\","), collapse = "")
          levelstring = substr(levelstring, 1, nchar(levelstring) - 1)
          finalstring = c(finalstring, "c(", levelstring, ")")
        }
      }
    }
    finalstring
  })

  regularmodelstring = reactive({
    tryCatch({
      if (any(unlist(strsplit(as.character(as.formula(input$model)[2]), "\\s\\+\\s|\\s\\*\\s|\\:")) == ".")) {
        dotreplace = paste0("(", paste0(names(candidatesetall()), collapse = " + "), ")")
        additionterms = unlist(strsplit(as.character(as.formula(input$model)[2]), "\\s\\+\\s"))
        multiplyterms = unlist(lapply(lapply(strsplit(additionterms, split = "\\s\\*\\s"), gsub, pattern = "^\\.$", replacement = dotreplace), paste0, collapse = " * "))
        interactionterms = unlist(lapply(lapply(strsplit(multiplyterms, split = "\\:"), gsub, pattern = "^\\.$", replacement = dotreplace), paste0, collapse = ":"))
        stringmodel = paste0("~", paste(interactionterms, collapse = " + "), sep = "")
      }
      paste0(as.character(as.formula(stringmodel)), collapse = "")
    }, error = function(e) {paste0(input$model, collapse = "")}
    )
  })

  modelwithblocks = reactive({
    if (isblockingtext()) {
      basemodel = gsub(pattern = "~", replacement = "", x = regularmodelstring(), fixed = TRUE)
      blockingmodelterms = "~ (1|Block1) + "
      paste0(blockingmodelterms, basemodel)
    }
  })

  contraststring = reactive({
    factorcat = c(input$factortype1, input$factortype2, input$factortype3, input$factortype4, input$factortype5, input$factortype6) == "cat"
    namecat = c(input$factorname1, input$factorname2, input$factorname3, input$factorname4, input$factorname5, input$factorname6)[factorcat]
    contrasttemp = "list("
    for (i in 1:length(namecat)) {
      if (i != length(namecat)) {
        contrasttemp = paste0(contrasttemp, namecat[i], " = ", "contr.sum, ")
      } else {
        contrasttemp = paste0(contrasttemp, namecat[i], " = ", "contr.sum)")
      }
    }
    contrasttemp
  })

  anyfactors = reactive({
    any(c(input$factortype1, input$factortype2, input$factortype3, input$factortype4, input$factortype5, input$factortype6) == "cat")
  })

  code = reactive({
    blocking = any("htc" %in% c(input$blockdepth1, input$blockdepth2, input$blockdepth3, input$blockdepth4, input$blockdepth5, input$blockdepth6)[1:input$numberfactors])

    first = paste0(c("<br><pre>",
                     "<code style=\"color:#468449\"># This is the R code used to generate these results in skpr.</code><br>",
                     "<code style=\"color:#468449\"># Copy this into an R script and rerun to reproduce these results.</code><br><br>",
                     "library(skpr)<br><br>",
                     ifelse(input$setseed,
                            paste0("<code style=\"color:#468449\">#Setting random number generator seed:</code><br>",
                                   "set.seed(", input$seed, ")<br><br>"), "<code style=\"color:#468449\">#Consider setting a seed to make this script fully reproducible.<br>#Go to Advanced->Set Random Number Generator Seed, click <br>#the checkbox, and set Random Seed to any whole number.</code><br><br>"),
                     "<code style=\"color:#468449\"># Generating candidate set:</code><br>",
                     "candidateset = expand.grid(", inputstring(), ")<br><br>",
                     ifelse(blocking,
                            paste0(c("<code style=\"color:#468449\"># Generating design for hard-to-change factors:</code> <br>",
                                     "design_htc = gen_design(candidateset = candidateset, <br>", rep("&nbsp;", 24),
                                     "model = ", blockmodel(), ", <br>", rep("&nbsp;", 24),
                                     "trials = ", as.character(input$numberblocks), ")<br><br>"), collapse = ""), ""),
                     "<code style=\"color:#468449\"># Generating design:</code><br>",
                     "design = gen_design(candidateset = candidateset, <br>", rep("&nbsp;", 20),
                     "model = ", regularmodelstring(), ", <br>", rep("&nbsp;", 20),
                     "trials = ", as.character(input$trials)
    ), collapse = "")
    if (blocking) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "splitplotdesign = design_htc"), collapse = "")
    }
    if (input$optimality != "D") {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "optimality = \"", input$optimality, "\""), collapse = "")
    }
    if (input$repeats != 20) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "repeats = ", input$repeats), collapse = "")
    }
    if (input$varianceratio != 1) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "varianceratio = ", input$varianceratio), collapse = "")
    }
    if (input$aliaspower != 2) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "aliaspower = ", input$aliaspower), collapse = "")
    }
    if (input$mindopt != 0.8) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "minDopt = ", input$mindopt), collapse = "")
    }
    if (as.logical(input$parallel)) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "parallel = TRUE"), collapse = "")
    }
    if (isblockingtext()) {
      first = paste(c(first, ", <br>", rep("&nbsp;", 20),
                      "add_blocking_columns = ", ifelse(input$splitanalyzable, "TRUE", "FALSE")), collapse = "")
    }
    first = paste0(c(first, ")<br><br>"), collapse = "")
    if (input$evaltype == "lm") {
      first = paste0(c(first,
                       "<code style=\"color:#468449\"># Evaluating Design:</code><br>",
                       "eval_design(design = design, <br>", rep("&nbsp;", 12),
                       "model = ", regularmodelstring(), ", <br>", rep("&nbsp;", 12),
                       "alpha = ", input$alpha), collapse = "")
      if (isblockingtext()) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 12),
                        "blocking = TRUE"), collapse = "")
      }
      if (input$snr != 2) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 12),
                        "effectsize = ", input$snr), collapse = "")
      }
      if (input$varianceratio != 1) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 12),
                        "varianceratios = ", input$varianceratio), collapse = "")
      }
      if (as.logical(input$conservative)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 12),
                        "conservative = TRUE"), collapse = "")
      }
      if (as.logical(input$detailedoutput)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 12),
                        "detailedoutput = TRUE"), collapse = "")
      }
      first = paste0(c(first, ")<br><br>"), collapse = "")
      first = paste0(first,
                     "<code style=\"color:#468449\">## How to analyze this experiment when the data have been collected:</code><br>",
                     "<code style=\"color:#468449\">## (to run, remove one # from this section) </code><br>",
                     "<code style=\"color:#468449\">## First, assign the results to a column in the data frame. Each </code><br>",
                     "<code style=\"color:#468449\">## entry in the vector corresponds to the result from that run in the design. <br><br></code>",
                     "<code style=\"color:#468449\">#design$Y = results <br><br></code>",
                     ifelse(!isblockingtext(),
                            "<code style=\"color:#468449\">## Now analyze the linear model with lm:</code><br><br>",
                            "<code style=\"color:#468449\">## Now analyze the blocked linear model with lmer (from the lme4 package) and lmerTest:<br><br></code>"),
                     ifelse(!isblockingtext(),
                            paste0("<code style=\"color:#468449\">#lm(formula = Y ", regularmodelstring(),
                                   ", data = design",
                                   ifelse(anyfactors(),
                                          paste0(", </code><br><code style=\"color:#468449\">#   ", "contrasts = ", contraststring(), ")</code>"),
                                          ")<br><br></code>")),
                            paste0(ifelse(input$splitanalyzable, "", "<code style=\"color:#468449\">## Note: Argument add_blocking_columns needs to be active in last gen_design call in order<br>## to analyze data taking into account the split-plot structure. The code below assumes that is true. <br><br></code>"),
                                   "<code style=\"color:#468449\">#library(lmerTest)<br>#lme4::lmer(formula = Y ",
                                   modelwithblocks(),
                                   ", data = design",
                                   ifelse(anyfactors(), paste0(", <br>#          ", "contrasts = ", contraststring(), "))<br><br>"), "))<br><br></code>"))))
    }
    if (input$evaltype == "glm") {
      first = paste0(c(first,
                       "<code style=\"color:#468449\"># Evaluating (Monte Carlo) Design:</code><br>",
                       "eval_design_mc(design = design, <br>", rep("&nbsp;", 15),
                       "model = ", regularmodelstring(), ", <br>", rep("&nbsp;", 15),
                       "alpha = ", input$alpha), collapse = "")
      if (isblockingtext()) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "blocking = TRUE"), collapse = "")
      }
      if (input$nsim != 1000) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "nsim = ", input$nsim), collapse = "")
      }
      if (input$glmfamily != "gaussian") {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "glmfamily = \"", input$glmfamily, "\""), collapse = "")
      }
      if (length(effectsize()) == 1) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "effectsize = ", effectsize()), collapse = "")
      } else {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "effectsize = c(", effectsize()[1], ", ", effectsize()[2], ")"), collapse = "")
      }
      if (!is.null(input$varianceratios)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "varianceratios = ", input$varianceratio), collapse = "")
      }
      if (as.logical(input$parallel_eval_glm)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "parallel = TRUE"), collapse = "")
      }
      if (as.logical(input$detailedoutput)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 15),
                        "detailedoutput = TRUE"), collapse = "")
      }
      first = paste0(c(first, ")<br><br>"), collapse = "")
      first = paste0(first,
                     "<code style=\"color:#468449\">## How to analyze this experiment when the data have been collected:</code><br>",
                     "<code style=\"color:#468449\">## (to run, remove one # from this section) </code><br>",
                     "<code style=\"color:#468449\">## First, assign the results to a column in the data frame. Each </code><br>",
                     "<code style=\"color:#468449\">## entry in the vector corresponds to the result from that run in the design. <br><br></code>",
                     "<code style=\"color:#468449\">#design$Y = results <br><br></code>",
                     ifelse(!isblockingtext(),
                            "<code style=\"color:#468449\">## Now analyze the generalized linear model with glm:</code><br><br>",
                            "<code style=\"color:#468449\">## Now analyze the blocked generalized linear model with glmer (from the lme4 package):<br><br></code>"),
                     ifelse(!isblockingtext(),
                            paste0("<code style=\"color:#468449\">#glm(formula = Y ", regularmodelstring(),
                                   ", data = design",
                                   ", <br>#   family = ", ifelse(input$glmfamily == "exponential", "Gamma(link=\"log\")", paste0("\"", input$glmfamily, "\"")),
                                   ifelse(anyfactors(),
                                          paste0(", </code><br><code style=\"color:#468449\">#   ", "contrasts = ", contraststring(), ")</code>"),
                                          ")<br><br></code>")),
                            paste0(ifelse(input$splitanalyzable, "", "<code style=\"color:#468449\">## Note: Argument add_blocking_columns needs to be active in last gen_design call in order<br>## to analyze data taking into account the split-plot structure. The code below assumes that is true. <br><br></code>"),
                                   "<code style=\"color:#468449\">#lme4::glmer(formula = Y ",
                                   modelwithblocks(),
                                   ", data = design",
                                   ", <br>#          family = ", ifelse(input$glmfamily == "exponential", "Gamma(link=\"log\")", paste0("\"", input$glmfamily, "\"")),
                                   ifelse(anyfactors(), paste0(", <br>#          ", "contrasts = ", contraststring(), ")"), ")</code>"))))
    }
    if (input$evaltype == "surv") {
      first = paste0(c(first,
                       "<code style=\"color:#468449\"># Evaluating (Monte Carlo Survival) Design:</code><br>",
                       "eval_design_survival_mc(design = design, <br>", rep("&nbsp;", 24),
                       "model = ", as.character(as.formula(input$model)), ", <br>", rep("&nbsp;", 24),
                       "alpha = ", input$alpha), collapse = "")
      if (input$nsim_surv != 1000) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "nsim = ", input$nsim_surv), collapse = "")
      }
      if (input$distribution != "gaussian") {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "distribution = \"", input$distribution, "\""), collapse = "")
      }
      if (length(effectsize()) == 1) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "effectsize = ", effectsize()), collapse = "")
      } else {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "effectsize = c(", effectsize()[1], ", ", effectsize()[2], ")"), collapse = "")
      }
      if (!is.na(input$censorpoint)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "censorpoint = ", input$censorpoint), collapse = "")
      }
      if (input$censortype != "right") {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "censortype = \"", input$censortype, "\""), collapse = "")
      }
      if (as.logical(input$parallel_eval_surv)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "parallel = TRUE"), collapse = "")
      }
      if (as.logical(input$detailedoutput)) {
        first = paste(c(first, ", <br>", rep("&nbsp;", 24),
                        "detailedoutput = TRUE"), collapse = "")
      }
      first = paste0(c(first, ")<br><br>"), collapse = "")
      first = paste0(first,
                     "<code style=\"color:#468449\">## How to analyze this experiment when the data have been collected:</code><br>",
                     "<code style=\"color:#468449\">## (to run, remove one # from this section) </code><br>",
                     "<code style=\"color:#468449\">## This is a survival model, so first create a Surv object that marks the censored runs.</code><br>",
                     "<code style=\"color:#468449\">## Each entry in the results vector corresponds to the result from that run in the design.</code><br>",
                     "<code style=\"color:#468449\">## Here, the raw censored responses are given as NA. We replace those values with the</code><br>",
                     "<code style=\"color:#468449\">## censor point and use the event argument to mark them as censored. </code><br>",
                     "<code style=\"color:#468449\">## (Surv argument \"event\" is TRUE when \"results\" is not censored, and FALSE when censored).<br><br></code>",
                     "<code style=\"color:#468449\">#notcensored = !is.na(results) </code><br>",
                     ifelse(!is.na(input$censorpoint), paste0("<code style=\"color:#468449\">#results[is.na(results)] = ", input$censorpoint, "</code><br>"),
                            ""),
                     "<code style=\"color:#468449\">#design$Y = Surv(time=results, event = notcensored, type = \"", input$censortype, "\") <br><br></code>",
                     "<code style=\"color:#468449\">## Now analyze the linear model with survreg (from the survival package):</code><br><br>",
                     paste0("<code style=\"color:#468449\">#survival::survreg(formula = Y ", regularmodelstring(),
                            ", data = design, dist = \"", input$distribution, "\")"))
    }

    first = paste0(c(first, "</pre></code>"), collapse = "")
    first
  })

  isblocking = reactive({
    input$submitbutton
    any("htc" %in% c(isolate(input$blockdepth1),
                     isolate(input$blockdepth2),
                     isolate(input$blockdepth3),
                     isolate(input$blockdepth4),
                     isolate(input$blockdepth5),
                     isolate(input$blockdepth6))[1:isolate(input$numberfactors)])
  })
  isblockingtext = reactive({
    any("htc" %in% c(input$blockdepth1,
                     input$blockdepth2,
                     input$blockdepth3,
                     input$blockdepth4,
                     input$blockdepth5,
                     input$blockdepth6)[1:(input$numberfactors)])
  })

  blockmodel = reactive({
    if (input$model == "~.") {
      as.formula(paste0("~", paste(names(inputlist_htctext()), collapse = " + ")))
    } else {
      names = names(inputlist())
      modelsplit = attr(terms.formula(as.formula(input$model), data = candidatesetall()), "term.labels")
      regularmodel = rep(FALSE, length(modelsplit))
      for (term in names) {
        regex = paste0("(\\b", term, "\\b)|(\\b", term, ":)|(:", term, "\\b)|(\\b", term, "\\s\\*)|(\\*\\s", term, "\\b)|(:", term, ":)")
        regularmodel = regularmodel | grepl(regex, modelsplit, perl = TRUE)
      }
      paste0("~", paste(modelsplit[!regularmodel], collapse = " + "))
    }
  })

  optimality = reactive({
    input$submitbutton
    if (isolate(input$numberfactors) == 1 && isolate(input$optimality) == "Alias") {
      showNotification("Alias-optimal design selected with only one factor: Switching to D-optimal.", type = "warning", duration = 10)
      updateSelectInput(session, "optimality", choices = c("D", "I", "A", "Alias", "G", "E", "T"), selected = "D")
      "D"
    } else {
      isolate(input$optimality)
    }
  })

  effectsize = reactive({
    if (input$evaltype == "lm") {
      return(input$snr)
    }
    if (input$evaltype == "glm") {
      if (input$glmfamily == "gaussian") {
        return(input$snr)
      }
      if (input$glmfamily == "binomial") {
        return(c(input$binomialprobs[1], input$binomialprobs[2]))
      }
      if (input$glmfamily == "poisson") {
        return( (c(input$poislow, input$poishigh)))
      }
      if (input$glmfamily == "exponential") {
        return(c(input$explow, input$exphigh))
      }
    }
    if (input$evaltype == "surv") {
      if (input$distribution == "gaussian") {
        return(input$snr)
      }
      if (input$distribution == "lognormal") {
        return(input$snr)
      }
      if (input$distribution == "exponential") {
        return(c(input$explow, input$exphigh))
      }
    }
  })

  runmatrix = reactive({
    input$submitbutton
    shinyjs::disable("submitbutton")
    shinyjs::disable("evalbutton")
    if (isolate(input$setseed)) {
      set.seed(isolate(input$seed))
    }
    tryCatch({
      if (!isblocking()) {
        withProgress(message = "Generating design:", value = 0, min = 0, max = 1, expr = {
          design = gen_design(candidateset = isolate(expand.grid(candidatesetall())),
                              model = isolate(as.formula(input$model)),
                              trials = isolate(input$trials),
                              optimality = isolate(optimality()),
                              repeats = isolate(input$repeats),
                              aliaspower = isolate(input$aliaspower),
                              minDopt = isolate(input$mindopt),
                              parallel = isolate(as.logical(input$parallel)),
                              advancedoptions = list(GUI = TRUE, progressBarUpdater = inc_progress_session))
        })
      } else {
        withProgress(message = "Generating whole-plots:", value = 0, min = 0, max = 1, expr = {
          spd = gen_design(candidateset = isolate(expand.grid(candidatesetall())),
                           model = isolate(as.formula(blockmodel())),
                           trials = isolate(input$numberblocks),
                           optimality = ifelse(toupper(isolate(optimality())) == "ALIAS" &&
                                                 length(isolate(inputlist_htc())) == 1, "D", isolate(optimality())),
                           repeats = isolate(input$repeats),
                           varianceratio = isolate(input$varianceratio),
                           aliaspower = isolate(input$aliaspower),
                           minDopt = isolate(input$mindopt),
                           parallel = isolate(as.logical(input$parallel)),
                           advancedoptions = list(GUI = TRUE, progressBarUpdater = inc_progress_session))
        })
        withProgress(message = "Generating full design:", value = 0, min = 0, max = 1, expr = {
          design = gen_design(candidateset = isolate(expand.grid(candidatesetall())),
                              model = isolate(as.formula(input$model)),
                              trials = isolate(input$trials),
                              splitplotdesign = spd,
                              optimality = isolate(optimality()),
                              repeats = isolate(input$repeats),
                              varianceratio = isolate(input$varianceratio),
                              aliaspower = isolate(input$aliaspower),
                              minDopt = isolate(input$mindopt),
                              parallel = isolate(as.logical(input$parallel)),
                              add_blocking_columns = isolate(input$splitanalyzable),
                              advancedoptions = list(GUI = TRUE, progressBarUpdater = inc_progress_session))
        })
      }
    }, finally = {
      shinyjs::enable("evalbutton")
      shinyjs::enable("submitbutton")
    })
    return(design)
  })

  evaluationtype = reactive({
    input$evalbutton
    isolate(input$evaltype)
  })

  format_table = function(powerval, display_table, alpha, nsim, colorblind) {
    color_bad = "red"
    color_maybe = "yellow"
    if(colorblind) {
      color_bad = "purple"
      color_maybe = "orange"
    }
    display_table = display_table %>%
      data_color(columns = "power",
                 palette = scales::col_numeric(palette =colorRampPalette(c("white", "darkgreen"))(100),
                                              domain =c(0,1)),
                 alpha = 0.3,
                 autocolor_text = FALSE) %>%
      tab_options(table.width = pct(100))
    if(any(powerval$power <= alpha + 1/sqrt(nsim) &
           powerval$power >= alpha - 1/sqrt(nsim))) {
      display_table = display_table %>%
        tab_style(
          style = list(
            cell_fill(color = color_maybe,alpha=0.3)
          ),
          locations = cells_body(
            columns = "power",
            rows = power <= alpha + 1/sqrt(nsim))
        ) %>%
        tab_source_note(
          source_note = sprintf("Note: Power values marked in %s are within the simulation uncertainty for user-specified Type-I error (increase the number of simulations)",
                                color_maybe)
        )
    }
    if(any(powerval$power < alpha - 1/sqrt(nsim))) {
      display_table = display_table %>%
        tab_style(
          style = list(
            cell_fill(color = color_bad,alpha=0.3)
          ),
          locations = cells_body(
            columns = "power",
            rows = (power < alpha - 1/sqrt(nsim)))
        ) %>%
        tab_source_note(
          source_note = sprintf("Note: Power values marked in %s fall below the user-specified Type-I error (%0.2f)",
                                color_bad, alpha)
        )
    }
    return(display_table)
  }

  powerresults = reactive({
    input$evalbutton
    if (evaluationtype() == "lm") {
      powerval = eval_design(design = isolate(runmatrix()),
                             model = as.formula(isolate(input$model)),
                             alpha = isolate(input$alpha),
                             blocking = isblocking(),
                             effectsize = isolate(effectsize()),
                             conservative = isolate(input$conservative),
                             detailedoutput = isolate(input$detailedoutput))
      powerval
    }
  })
  powerresultsglm = reactive({
    input$evalbutton
    if (isolate(evaluationtype()) == "glm") {
      if (isolate(input$setseed)) {
        set.seed(isolate(input$seed))
      }
      withProgress(message = ifelse(isolate(isblocking()), "Simulating (with REML):", "Simulating:"), value = 0, min = 0, max = 1, expr = {
        powerval = suppressWarnings(eval_design_mc(design = isolate(runmatrix()),
                                                   model = isolate(as.formula(input$model)),
                                                   alpha = isolate(input$alpha),
                                                   blocking = isolate(isblocking()),
                                                   nsim = isolate(input$nsim),
                                                   varianceratios = isolate(input$varianceratio),
                                                   glmfamily = isolate(input$glmfamily),
                                                   effectsize = isolate(effectsize()),
                                                   parallel = isolate(input$parallel_eval_glm),
                                                   detailedoutput = isolate(input$detailedoutput),
                                                   advancedoptions = list(GUI = TRUE, progressBarUpdater = inc_progress_session)))
      })
      powerval
    }
  })
  powerresultssurv = reactive({
    input$evalbutton
    if (isolate(evaluationtype()) == "surv") {
      if (isolate(input$setseed)) {
        set.seed(isolate(input$seed))
      }
      if (isolate(isblocking())) {
        print("Hard-to-change factors are not supported for survival designs. Evaluating design with no blocking.")
      }
      withProgress(message = "Simulating:", value = 0, min = 0, max = 1, expr = {
        powerval = suppressWarnings(eval_design_survival_mc(design = isolate(runmatrix()),
                                                            model = isolate(as.formula(input$model)),
                                                            alpha = isolate(input$alpha),
                                                            nsim = isolate(input$nsim_surv),
                                                            censorpoint = isolate(input$censorpoint),
                                                            censortype = isolate(input$censortype),
                                                            distribution = isolate(input$distribution),
                                                            effectsize = isolate(effectsize()),
                                                            detailedoutput = isolate(input$detailedoutput),
                                                            advancedoptions = list(GUI = TRUE, progressBarUpdater = inc_progress_session)))
        powerval
      })
    }
  })

  pal_option = function(n) {
    if(input$colorchoice == "A") {
      viridis::magma(n)
    } else if(input$colorchoice == "B") {
      viridis::inferno(n)
    } else if(input$colorchoice == "C") {
      viridis::plasma(n)
    } else if(input$colorchoice == "D") {
      viridis::viridis(n)
    } else {
      "white"
    }
  }

  style_matrix = function(runmat, order_vals = FALSE, alpha = 0.3, trials, optimality) {
    . = NULL
    if(order_vals) {
      new_runmat = runmat[do.call(order, runmat),, drop=FALSE ]
      rownames(new_runmat) = 1:nrow(new_runmat)

      display_rm = gt(new_runmat[do.call(order, new_runmat),, drop=FALSE ],
                      rownames_to_stub = TRUE) %>%
        tab_stubhead("Run") %>%
        tab_options(data_row.padding = px(10))  %>%
        tab_spanner(
          label = "Factors",
          columns = colnames(.)
        ) %>% tab_header(
          title = "Design",
          subtitle = sprintf("%d-run %s-optimal design",
                             trials,
                             optimality)
        )
    } else {
      display_rm = gt(runmat,rownames_to_stub = TRUE) %>%
        tab_stubhead("Run") %>%
        tab_options(data_row.padding = px(10)) %>%
        tab_spanner(
          label = "Factors",
          columns = colnames(.)
        ) %>% tab_header(
          title = "Design",
          subtitle = sprintf("%d-run %s-optimal design",
                             trials,
                             optimality)
        )
    }
    cols_rm = colnames(runmat)
    for(i in seq_len(length(cols_rm))) {
      if(is.numeric(runmat[,i])) {
        display_rm = display_rm %>%
          data_color(
            columns = cols_rm[i],
            palette = pal_option(100),
            alpha = alpha,
            autocolor_text = FALSE)
      } else {
        display_rm = display_rm %>%
          data_color(
            columns = cols_rm[i],
            palette = pal_option(length(unique(runmat[,i]))),
            alpha = alpha,
            autocolor_text = FALSE)
      }
    }
    display_rm
  }

  output$runmatrix = gt::render_gt({
    style_matrix(runmatrix(), order_vals = input$orderdesign,  trials = isolate(input$trials), optimality = isolate(input$optimality))
  }, align = "left")

  output$powerresults = gt::render_gt( {
    req(powerresults())
    format_table(powerresults(),gt(powerresults()), isolate(input$alpha),isolate(input$nsim),isolate(input$colorblind))
  }, align = "left")

  output$powerresultsglm = gt::render_gt( {
    req(powerresultsglm())

    format_table(powerresultsglm(),gt(powerresultsglm()), isolate(input$alpha),isolate(input$nsim),isolate(input$colorblind))
  }, align = "left")

  output$powerresultssurv = gt::render_gt({
    req(powerresultssurv())

    format_table(powerresultssurv(),gt(powerresultssurv()), isolate(input$alpha),isolate(input$nsim_surv),isolate(input$colorblind))
  }, align = "left")

  output$aliasplot = renderPlot({
    input$submitbutton
    tryCatch({
      plot_correlations(isolate(runmatrix()))
    }, error = function(e) {
    })
  })

  output$fdsplot = renderPlot({
    input$submitbutton
    plot_fds(isolate(runmatrix()))
  })

  output$code = renderUI({
    HTML(code())
  })

  output$dopt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "D"))
  })
  output$aopt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "A"))
  })
  output$iopt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "I"))
  })
  output$eopt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "E"))
  })
  output$gopt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "G"))
  })
  output$topt = renderText({
    input$submitbutton
    isolate(attr(runmatrix(), "T"))
  })
  output$optimalsearch = renderPlot({
    input$submitbutton
    if (isolate(optimality()) %in% c("D", "G", "A")) {
      if(attr(runmatrix(), "blocking") || attr(runmatrix(), "splitplot")) {
        max_y_val = max(attr(runmatrix(), "optimalsearchvalues"),na.rm=TRUE)
        statement = "Optimality Value (higher is better)"
      }  else {
        max_y_val = 100
        statement = "Efficiency (higher is better)"
      }
      isolate(plot(attr(runmatrix(), "optimalsearchvalues"), xlab = "Search Iteration", ylab = paste(optimality(), statement),
                   type = "p", col = "red", pch = 16, ylim = c(0, max_y_val)))
      isolate(points(x = attr(runmatrix(), "best"), y = attr(runmatrix(), "optimalsearchvalues")[attr(runmatrix(), "best")],
                     type = "p", col = "green", pch = 16, cex = 2, ylim = c(0, max_y_val)))
    } else {
      if (isolate(optimality()) == "I") {
        isolate(plot(attr(runmatrix(), "optimalsearchvalues"), xlab = "Search Iteration", ylab = "Average Prediction Variance (lower is better)", type = "p", col = "red", pch = 16))
      } else {
        isolate(plot(attr(runmatrix(), "optimalsearchvalues"), xlab = "Search Iteration", ylab = paste(optimality(), "Criteria Value (higher is better)"), type = "p", col = "red", pch = 16))
      }
      isolate(points(x = attr(runmatrix(), "best"), y = attr(runmatrix(), "optimalsearchvalues")[attr(runmatrix(), "best")], type = "p", col = "green", pch = 16, cex = 2))
    }
  })
  output$simulatedpvalues = renderPlot({
    updateval = c(powerresultsglm(),powerresultssurv())
    if(isolate(evaluationtype() == "glm")) {
      pvalrows = isolate(floor(ncol(attr(powerresultsglm(), "pvals")) / 3) + 1)
      if (!is.null(attr(powerresultsglm(), "pvals"))) {
        par(mfrow = c(pvalrows, 3))
        for (col in 1:isolate(ncol(attr(powerresultsglm(), "pvals")))) {
          isolate(hist(attr(powerresultsglm(), "pvals")[, col], breaks = seq(0, 1, 0.05),
                       main = colnames(attr(powerresultsglm(), "pvals"))[col],
                       xlim = c(0, 1), xlab = "p values", ylab = "Count", col = "red", pch = 16))
        }
      }
    }
    if(isolate(evaluationtype() == "surv")) {
      pvalrows = isolate(floor(ncol(attr(powerresultssurv(), "pvals")) / 3) + 1)
      if (!is.null(attr(powerresultssurv(), "pvals"))) {
        par(mfrow = c(pvalrows, 3))
        for (col in 1:isolate(ncol(attr(powerresultssurv(), "pvals")))) {
          isolate(hist(attr(powerresultssurv(), "pvals")[, col], breaks = seq(0, 1, 0.05),
                       main = colnames(attr(powerresultssurv(), "pvals"))[col],
                       xlim = c(0, 1), xlab = "p values", ylab = "Count", col = "red", pch = 16))
        }
      }
    }
  })
  output$parameterestimates = renderPlot({
    input$evalbutton
    if (!is.null(attr(powerresultsglm(), "estimates"))) {
      ests = apply(attr(powerresultsglm(), "estimates"), 2, quantile, c(0.05, 0.5, 0.95))
      truth = attr(powerresultsglm(), "anticoef")
      if (isolate(input$glmfamily) == "binomial") {
        ests = exp(ests) / (1 + exp(ests))
        truth = exp(truth) / (1 + exp(truth))
      }
      if (isolate(input$glmfamily) == "poisson") {
        ests = exp(ests)
        truth = exp(truth)
      }
      if (isolate(input$glmfamily) == "exponential") {
        ests = exp(ests)
        truth = exp(truth)
      }
      par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
      plot(x = 1:length(colnames(ests)), y = ests[2, ],
           xaxt = "n",
           xlab = "Parameters",
           ylab = ifelse(isolate(input$glmfamily) == "binomial", "Parameter Estimates (Probability)", "Parameter Estimates"),
           ylim = ifelse(rep(isolate(input$glmfamily) == "binomial", 2), c(0, 1), c(min(as.vector(ests)), max(as.vector(ests)))),
           xlim = c(0.5, length(colnames(ests)) + 0.5),
           type = "p", pch = 16, col = "red", cex = 1)
      axis(1, at = 1:length(colnames(ests)), labels = colnames(ests), las = 2)
      legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
      par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
      grid(nx = NA, ny = NULL)
      arrows(x0 = 1:length(colnames(ests)), y0 = ests[1, ], x1 = 1:length(colnames(ests)), y1 = ests[3, ], length = 0.05, angle = 90, code = 3)
      points(x = 1:length(colnames(ests)), y = truth, pch = 16, col = "blue", cex = 1)
      title("Simulated Parameter Estimates (5%-95% Confidence Intervals)")
    }
  })

  output$parameterestimatessurv = renderPlot({
    input$evalbutton
    if (!is.null(attr(powerresultssurv(), "estimates"))) {
      ests = apply(attr(powerresultssurv(), "estimates"), 2, quantile, c(0.05, 0.5, 0.95))
      truth = attr(powerresultssurv(), "anticoef")
      if (isolate(input$distibution) == "exponential") {
        ests = exp(ests)
        truth = exp(truth)
      }
      par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
      plot(x = 1:length(colnames(ests)), y = ests[2, ],
           xaxt = "n",
           xlab = "Parameters",
           ylab = ifelse(isolate(input$glmfamily) == "binomial", "Parameter Estimates (Probability)", "Parameter Estimates"),
           ylim = ifelse(rep(isolate(input$glmfamily) == "binomial", 2), c(0, 1), c(min(as.vector(ests)), max(as.vector(ests)))),
           xlim = c(0.5, length(colnames(ests)) + 0.5),
           type = "p", pch = 16, col = "red", cex = 1)
      axis(1, at = 1:length(colnames(ests)), labels = colnames(ests), las = 2)
      legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
      par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
      grid(nx = NA, ny = NULL)
      arrows(x0 = 1:length(colnames(ests)), y0 = ests[1, ], x1 = 1:length(colnames(ests)), y1 = ests[3, ], length = 0.05, angle = 90, code = 3)
      points(x = 1:length(colnames(ests)), y = truth, pch = 16, col = "blue", cex = 1)
      title("Simulated Parameter Estimates (5%-95% Confidence Intervals)")
    }
  })

  output$responsehistogram = renderPlot({
    input$evalbutton
    if (!is.null(attr(powerresultsglm(), "estimates"))) {
      responses = as.vector(attr(powerresultsglm(), "estimates") %*% t(attr(powerresultsglm(), "modelmatrix")))
      trueresponses = as.vector(attr(powerresultsglm(), "anticoef") %*% t(attr(powerresultsglm(), "modelmatrix")))
      widths = hist(trueresponses, plot = FALSE)$counts
      widths = widths[widths != 0]
      widths = sqrt(widths)
      uniquevalues = length(table(responses))
      breakvalues = ifelse(uniquevalues < isolate(input$nsim) * isolate(input$trials) / 10, uniquevalues, isolate(input$nsim) * isolate(input$trials) / 10)
      if (isolate(input$glmfamily) == "binomial") {
        responses = exp(responses) / (1 + exp(responses))
        trueresponses = exp(trueresponses) / (1 + exp(trueresponses))
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response (Probability)", main = "Distribution of Simulated Response Estimates", xlim = c(0, 1))
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses Estimates", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)
      }
      if (isolate(input$glmfamily) == "poisson") {
        responses = exp(responses)
        trueresponses = exp(trueresponses)
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response", main = "Distribution of Simulated Response Estimates", xlim = c(ifelse(is.na(input$estimatesxminglm), min(hist(responses, plot = FALSE)$breaks), input$estimatesxminglm), ifelse(is.na(input$estimatesxmaxglm), max(hist(responses, plot = FALSE)$breaks), input$estimatesxmaxglm)), col = "red", border = "red")
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses ", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)
      }
      if (isolate(input$glmfamily) == "exponential") {
        responses = exp(responses)
        trueresponses = exp(trueresponses)
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response", main = "Distribution of Simulated Response Estimates", xlim = c(ifelse(is.na(input$estimatesxminglm), min(hist(responses, plot = FALSE)$breaks), input$estimatesxminglm), ifelse(is.na(input$estimatesxmaxglm), max(hist(responses, plot = FALSE)$breaks), input$estimatesxmaxglm)), col = "red", border = "red")
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)
      }
      if (isolate(input$glmfamily) == "gaussian") {
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response", main = "Distribution of Simulated Response Estimates", xlim = c(ifelse(is.na(input$estimatesxminglm), min(hist(responses, plot = FALSE)$breaks), input$estimatesxminglm), ifelse(is.na(input$estimatesxmaxglm), max(hist(responses, plot = FALSE)$breaks), input$estimatesxmaxglm)), col = "red", border = "red")
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)

      }
    }
  })

  output$responsehistogramsurv = renderPlot({
    input$evalbutton
    if (!is.null(attr(powerresultssurv(), "estimates"))) {
      responses = as.vector(attr(powerresultssurv(), "estimates") %*% t(attr(powerresultssurv(), "modelmatrix")))
      trueresponses = as.vector(attr(powerresultssurv(), "anticoef") %*% t(attr(powerresultssurv(), "modelmatrix")))
      filtered_string = ""
      if(isolate(input$distribution) == "exponential") {
        #Filter out extreme values
        mad_trueresp = 10*max(exp(trueresponses))
        num_filtered = sum(exp(responses) > mad_trueresp)
        responses = responses[exp(responses) < mad_trueresp]
        trueresponses = trueresponses[exp(trueresponses) < mad_trueresp]
        filtered_string = sprintf(" (%g extreme outliers removed)",num_filtered)
      }
      widths = hist(trueresponses, plot = FALSE)$counts
      widths = widths[widths != 0]
      widths = sqrt(widths)
      uniquevalues = length(table(responses))
      breakvalues = ifelse(uniquevalues < isolate(input$nsim_surv) * isolate(input$trials) / 10, uniquevalues, isolate(input$nsim_surv) * isolate(input$trials) / 10)
      if (isolate(input$distribution) == "exponential") {
        responses = exp(responses)
        trueresponses = exp(trueresponses)
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response", main = sprintf("Distribution of Simulated Response Estimates%s", filtered_string), xlim = c(ifelse(is.na(input$estimatesxminsurv), min(hist(responses, plot = FALSE)$breaks), input$estimatesxminsurv), ifelse(is.na(input$estimatesxmaxsurv), max(hist(responses, plot = FALSE)$breaks), input$estimatesxmaxsurv)), col = "red", border = "red")
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)
      }
      if (isolate(input$distribution) %in% c("gaussian", "lognormal")) {
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
        hist(responses, breaks = breakvalues, xlab = "Response", main = "Distribution of Simulated Response Estimates (from survival analysis)", xlim = c(ifelse(is.na(input$estimatesxminsurv), min(hist(responses, plot = FALSE)$breaks), input$estimatesxminsurv), ifelse(is.na(input$estimatesxmaxsurv), max(hist(responses, plot = FALSE)$breaks), input$estimatesxmaxsurv)), col = "red", border = "red")
        legend("topright", inset = c(-0.2, 0), legend = c("Truth", "Simulated"), pch = c(16, 16), col = c("blue", "red"), title = "Estimates")
        par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = FALSE)
        grid(nx = NA, ny = NULL)
        hist(responses, breaks = breakvalues, add = TRUE, main = "Distribution of Simulated Responses", xlab = "Response", ylab = "Count", col = "red", border = "red")
        abline(v = unique(trueresponses)[order(unique(trueresponses))], col = adjustcolor("blue", alpha.f = 0.40), lwd = widths)
      }
    }
  })
  output$separationwarning = renderText({
    input$evalbutton
    likelyseparation = FALSE
    if (isolate(input$evaltype) == "glm" && isolate(input$glmfamily) == "binomial") {
      if (!is.null(attr(powerresultsglm(), "pvals"))) {
        pvalmat = attr(powerresultsglm(), "pvals")
        for (i in 2:ncol(pvalmat)) {
          pvalcount = hist(pvalmat[, i], breaks = seq(0, 1, 0.05), plot = FALSE)
          likelyseparation = likelyseparation || (all(pvalcount$count[20] > pvalcount$count[17:19]) && pvalcount$count[20] > isolate(input$nsim) / 15)
        }
      }
    }
    if (likelyseparation) {
      showNotification("Partial or complete separation likely detected in the binomial Monte Carlo simulation. Increase the number of runs in the design or decrease the number of model parameters to improve power.", type = "warning", duration = 10)
    }
  })
  observeEvent(input$tutorial,
               introjs(session,
                       options = list("showProgress" = "true",
                                      "showBullets" = "false"),
                       events = list(
                         "onchange" = I("
                                          if (this._currentStep==0) {
                                          $('a[data-value=\"Advanced\"]').removeClass('active');
                                          $('a[data-value=\"Power\"]').removeClass('active');
                                          $('a[data-value=\"Basic\"]').addClass('active');
                                          $('a[data-value=\"Basic\"]').trigger('click');
                                          }
                                          if (this._currentStep==5) {
                                          $('a[data-value=\"Power\"]').removeClass('active');
                                          $('a[data-value=\"Basic\"]').removeClass('active');
                                          $('a[data-value=\"Advanced\"]').addClass('active');
                                          $('a[data-value=\"Advanced\"]').trigger('click');
                                          }
                                          if (this._currentStep==13) {
                                          $('a[data-value=\"Advanced\"]').removeClass('active');
                                          $('a[data-value=\"Power\"]').addClass('active');
                                          $('a[data-value=\"Power\"]').trigger('click');
                                          }
                                          if (this._currentStep==17) {
                                          $('input[value=\"glm\"]').trigger('click');
                                          }
                                          if (this._currentStep==21) {
                                          $('input[value=\"surv\"]').trigger('click');
                                          }
                                          if (this._currentStep==24) {
                                          $('a[data-value=\"Design Evaluation\"]').removeClass('active');
                                          $('a[data-value=\"Generating Code\"]').removeClass('active');
                                          $('a[data-value=\"Design\"]').addClass('active');
                                          $('a[data-value=\"Design\"]').trigger('click');
                                          $('#evaltype').val('glm');
                                          Shiny.onInputChange('evaltype', 'glm');
                                          $('#numberfactors').val('3');
                                          Shiny.onInputChange('numberfactors', 3);
                                          $('#trials').val('12');
                                          Shiny.onInputChange('trials', 12);
                                          $('#submitbutton').trigger('click');
                                          $('#evalbutton').trigger('click');
                                          }
                                          if (this._currentStep==25) {
                                          $('#evalbutton').trigger('click');
                                          $('a[data-value=\"Design\"]').removeClass('active');
                                          $('a[data-value=\"Design Evaluation\"]').addClass('active');
                                          $('a[data-value=\"Design Evaluation\"]').trigger('click');
                                          }
                                          if (this._currentStep==31) {
                                          $('a[data-value=\"Design Evaluation\"]').removeClass('active');
                                          $('a[data-value=\"Generating Code\"]').addClass('active');
                                          $('a[data-value=\"Generating Code\"]').trigger('click');
                                          }"
                         ))
               ))
  outputOptions(output, "separationwarning", suspendWhenHidden = FALSE)
}