# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# DIF/FAIRNESS ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####

# ** DIF matching variable (DMV) presence control ####
DIF_matching_presence <- reactive({
  !(length(dataset$DIFmatching) == 1 & any(dataset$DIFmatching == "missing"))
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * TOTAL SCORES ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

# ** Matching variable selection ####
# total scores or uploaded DIF matching variable
DIF_total_matching <- reactive({
  switch(input$DIF_total_matching,
    "score" = total_score(),
    # "zscore" = scale(apply(as.data.frame(total_score()), 1, sum)),
    "uploaded" = unlist(DIFmatching())
    # "zuploaded" = scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  )
})

# ** Title update ####
# one reactive output cannot be displayed more than once for some reason
output$DIF_total_matching_title1 <-
  output$DIF_total_matching_title2 <-
  output$DIF_total_matching_title3 <-
  output$DIF_total_matching_title4 <-
  output$DIF_total_matching_title5 <-
  output$DIF_total_matching_title6 <-
  renderText({
    ifelse(input$DIF_total_matching == "uploaded",
      "uploaded variable",
      "total scores"
    )
  })

# ** Updating options for DIF matching variable ####
observe({
  if (DIF_matching_presence()) {
    updateSelectInput(
      session,
      "DIF_total_matching",
      choices = c(
        "Total score" = "score",
        # "Standardized total score" = "zscore",
        "Uploaded" = "uploaded"
        # "Standarized uploaded" = "zuploaded"
      ),
      selected = "score"
    )
  } else {
    updateSelectInput(
      session,
      "DIF_total_matching",
      choices = c(
        "Total score" = "score"
        # "Standardized total score" = "zscore"
      ),
      selected = "score"
    )
  }
})

# ** Tooltips  ####
# copy of Total scores from Scores
output$DIF_total_matching_tooltip_mean <- renderUI({
  span(
    class = "ttooltip",
    style = "color: #2286bf",
    "mean,",
    span(
      class = "ttooltiptext",
      "$$\\bar{X} = \\frac{1}{n}\\sum_{p = 1}^n X_p$$"
    )
  )
})
output$DIF_total_matching_tooltip_sd <- renderUI({
  span(
    class = "ttooltip",
    style = "color: #2286bf",
    "sample standard deviation",
    span(
      class = "ttooltiptext",
      "$$\\textrm{SD}(X) = \\sqrt{\\frac{1}{n - 1}\\sum_{p = 1}^n (X_p - \\bar{X})^2}$$"
    )
  )
})
output$DIF_total_matching_tooltip_skewness <- renderUI({
  span(
    class = "ttooltip",
    style = "color: #2286bf",
    "sample skewness,",
    span(
      class = "ttooltiptext",
      "$$\\frac{\\frac{1}{n} \\sum_{p = 1}^n (X_p - \\bar{X})^3}{\\left[\\frac{1}{n} \\sum_{p = 1}^n (X_p - \\bar{X})^2\\right]^{3/2}}$$"
    )
  )
})
output$DIF_total_matching_tooltip_kurtosis <- renderUI({
  span(
    class = "ttooltip",
    style = "color: #2286bf",
    "sample kurtosis.",
    span(
      class = "ttooltiptext",
      "$$\\frac{\\frac{1}{n} \\sum_{p = 1}^n (X_p - \\bar{X})^4}{\\left[\\frac{1}{n} \\sum_{p = 1}^n (X_p - \\bar{X})^2\\right]^2}$$"
    )
  )
})

# ** DIF matching variable with group variable ####
DIF_total_score_with_group <- reactive({
  tibble(
    total_matching = DIF_total_matching(),
    Group = factor(group(), levels = c(0, 1), labels = c("Reference (0)", "Focal (1)"))
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MAIN TAB ####

# ** Summary of total scores for groups ####
DIF_total_table <- reactive({
  d <- DIF_total_score_with_group()

  d |>
    group_by(Group) |>
    summarise(
      across(
        total_matching,
        list(
          n = length,
          `n<sub>c</sub>` = ~ sum(!is.na(.x)),
          Min = ~ min(.x, na.rm = TRUE),
          Max = ~ max(.x, na.rm = TRUE),
          Mean = ~ mean(.x, na.rm = TRUE),
          Median = ~ median(.x, na.rm = TRUE),
          SD = ~ sd(.x, na.rm = TRUE),
          Skewness = ~ ShinyItemAnalysis:::skewness(.x),
          Kurtosis = ~ ShinyItemAnalysis:::kurtosis(.x)
        ),
        .names = "{.fn}"
      )
    )
})

output$DIF_total_table <- renderTable(
  DIF_total_table(),
  digits = 2,
  sanitize.colnames.function = function(x) x
)

# ** Histogram of total scores ####
DIF_total_histogram <- reactive({
  xlab <- switch(input$DIF_total_matching,
    "score" = "Total score",
    "uploaded" = "Observed score"
  )

  d <- DIF_total_score_with_group()

  d |>
    ggplot(aes(x = total_matching, fill = Group, col = Group)) +
    # text aes to be picked by ggplotly and used as a tooltip
    geom_histogram(
      aes(
        text =
          paste0(
            "Count: ", after_stat(count), "\n",
            xlab, ": ", x, "\n",
            "Group: ", fill
          )
      ),
      binwidth = 1, position = "dodge2", alpha = 0.75
    ) +
    xlab(xlab) +
    ylab("Number of respondents") +
    scale_fill_manual(
      values = c(`Reference (0)` = "dodgerblue2", `Focal (1)` = "goldenrod2"),
      labels = c(`Reference (0)` = "Reference", `Focal (1)` = "Focal")
    ) +
    scale_colour_manual(
      values = c(`Reference (0)` = "dodgerblue2", `Focal (1)` = "goldenrod2"),
      labels = c(`Reference (0)` = "Reference", `Focal (1)` = "Focal")
    ) +
    theme_app()
})

output$DIF_total_histogram <- renderPlotly({
  DIF_total_histogram() |>
    ggplotly(tooltip = "text") |>
    plotly::config(displayModeBar = FALSE)
})

# ** Download histogram of total scores ####
output$DIF_total_histogram_download <- downloadHandler(
  filename = function() {
    "fig_HistogramGroups.png"
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_total_histogram() +
        theme(
          text = element_text(size = setting_figures$text_size),
          legend.position = "inside",
          legend.position.inside = c(0.15, 0.85)
        ),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** t-test to compare total scores ####
DIF_total_ttest <- reactive({
  sc0 <- DIF_total_matching()[group() == 0]
  sc1 <- DIF_total_matching()[group() == 1]

  validate(
    need(length(na.omit(sc0)) > 1L, "Not enough complete observations in the reference group."),
    need(length(na.omit(sc1)) > 1L, "Not enough complete observations in the focal group."),
    errorClass = "validation-error"
  )

  ttest <- t.test(sc0, sc1)

  tab <- c(
    paste0(
      sprintf("%.2f", mean(sc0, na.rm = TRUE) - mean(sc1, na.rm = TRUE)), " (",
      sprintf("%.2f", ttest$conf.int[1]), ", ",
      sprintf("%.2f", ttest$conf.int[2]), ")"
    ),
    sprintf("%.2f", ttest$statistic),
    sprintf("%.2f", ttest$parameter),
    ifelse(ttest$p.value < 0.001, "< 0.001", sprintf("%.3f", ttest$p.value))
  )
  tab <- t(as.data.frame(tab))
  colnames(tab) <- c("Diff. (CI)", "t-value", "df", "p-value")
  tab
})

output$DIF_total_ttest <- renderTable({
  DIF_total_ttest()
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * DELTA PLOT ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MODEL ####

# ** Delta plot calculations ####
DIF_DP <- reactive({
  data <- data.frame(binary(), group = group())
  switch(input$DIF_DP_threshold,
    "Fixed" = deltaPlot(data,
      group = "group",
      focal.name = 1,
      thr = 1.5,
      purify = input$DIF_DP_purification
      # purType = input$DIF_DP_purification_type
      # (only 1IPP is used, so no need to pass this arg)
    ),
    "Normal" = deltaPlot(data,
      group = "group",
      focal.name = 1,
      thr = "norm",
      purify = input$DIF_DP_purification,
      purType = input$DIF_DP_purification_type
    )
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MAIN TAB ####

# ** Delta plot ####
DIF_DP_plot <- reactive({
  dp <- DIF_DP()
  df <- data.frame(dp$Deltas)
  df$nam <- item_numbers()

  par <- dp$axis.par
  thr <- dp$thr

  if (length(par) > 2) {
    par <- par[length(par) / 2, ]
  }

  if (length(thr) > 1) {
    thr <- thr[length(thr)]
  }

  p <- ggplot(
    df,
    aes(x = X1, y = X2, label = nam)
  ) +
    geom_point() +
    geom_text(hjust = 0, nudge_x = 0.05, size = 6) +
    geom_abline(
      intercept = par[1], slope = par[2],
      linewidth = 1
    ) +
    geom_abline(
      intercept = par[1] + thr * sqrt(par[2]^2 + 1),
      slope = par[2],
      color = "red",
      linetype = "dashed",
      linewidth = 1
    ) +
    geom_abline(
      intercept = par[1] - thr * sqrt(par[2]^2 + 1),
      slope = par[2],
      color = "red",
      linetype = "dashed",
      linewidth = 1
    ) +
    labs(
      x = "Reference group",
      y = "Focal group"
    ) +
    xlim(min(dp$Deltas, na.rm = TRUE) - 0.5, max(dp$Deltas, na.rm = TRUE) + 0.5) +
    ylim(min(dp$Deltas, na.rm = TRUE) - 0.5, max(dp$Deltas, na.rm = TRUE) + 0.5) +
    theme_app()
  if (is.numeric(dp$DIFitems)) {
    df2 <- df[dp$DIFitems, ]
    p <- p + geom_point(
      data = df2,
      aes(x = X1, y = X2, label = nam),
      size = 8, color = "black", shape = 1
    )
  }
  p <- p + ggtitle("Delta plot")
  p
})

output$DIF_DP_plot <- renderPlotly({
  g <- DIF_DP_plot() |>
    style(textposition = "right")
  p <- ggplotly(g)

  p$x$data[[1]]$text <- gsub("X1", "Delta - reference", p$x$data[[1]]$text)
  p$x$data[[1]]$text <- gsub("X2", "Delta - focal", p$x$data[[1]]$text)
  p$x$data[[1]]$text <- gsub("nam:", "Item:", p$x$data[[1]]$text)
  p$x$data[[2]]$hovertext <- gsub("X1", "Delta - reference", p$x$data[[2]]$hovertext)
  p$x$data[[2]]$hovertext <- gsub("X2", "Delta - focal", p$x$data[[2]]$hovertext)
  p$x$data[[2]]$hovertext <- gsub("nam:", "Item:", p$x$data[[2]]$hovertext)

  for (i in 3:5) {
    p$x$data[[i]]$text <- gsub("intercept", "Intercept", p$x$data[[i]]$text)
    p$x$data[[i]]$text <- gsub("slope", "Slope", p$x$data[[i]]$text)
  }

  p$elementId <- NULL
  p |> plotly::config(displayModeBar = FALSE)
})

# ** Download delta plot ####
output$DIF_DP_plot_download <- downloadHandler(
  filename = function() {
    "fig_DIF_DP.png"
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_DP_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Summary table ####
# output table for the LAST iteration
DIF_DP_table <- reactive({
  res <- DIF_DP()
  ncR <- formatC(colSums(!is.na(binary()[group() == 0])), digits = 0, format = "f")
  ncF <- formatC(colSums(!is.na(binary()[group() == 1])), digits = 0, format = "f")
  # from print.deltaPlot:
  tab <- data.frame(ncR, ncF, round(cbind(res$Props, res$Deltas, res$Dist[, ncol(res$Dist)]), 4))
  symb <- symnum(abs(as.numeric(res$Dist[, ncol(res$Dist)])),
    c(0, abs(res$thr[length(res$thr)]), Inf),
    symbols = c("", "***")
  )
  tab$symb <- symb

  colnames(tab) <- c(
    "\\(\\mathit{n_{c, ref}}\\)", "\\(\\mathit{n_{c, foc}}\\)",
    "Prop. correct \\(\\mathit{_{ref}}\\)", "Prop. correct \\(\\mathit{_{foc}}\\)",
    "\\(\\mathit{\\Delta_{ref}}\\)", "\\(\\mathit{\\Delta_{foc}}\\)",
    "Distance", ""
  )
  rownames(tab) <- item_names()
  tab
})

output$DIF_DP_coef <- renderTable(
  {
    DIF_DP_table()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Download summary table ####
output$DIF_DP_table_download <- downloadHandler(
  filename = function() {
    "DIF_DP_statistics.csv"
  },
  content = function(file) {
    data <- DIF_DP_table()
    colnames(data) <- c(
      "n  (Ref.)", "n  (Foc.)",
      "Ref. prop. correct",
      "Foc. prop. correct",
      "Delta (Ref.)",
      "Delta (Foc.)",
      "Distance",
      ""
    )

    write.csv(data, file)
    write(paste(
      "Note:",
      DIF_DP_table_note()$puri,
      DIF_DP_table_note()$axes,
      DIF_DP_table_note()$thr,
      "Items detected as DIF are flagged with '***'",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Note about setting below summary table ####
DIF_DP_table_note <- reactive({
  res <- NULL

  model <- DIF_DP()
  thr <- model$thr
  thr <- thr[length(thr)]

  axis <- round(model$axis.par, 2)
  last <- nrow(axis)
  if (is.null(last)) {
    a <- axis[1]
    b <- axis[2]
  } else {
    a <- axis[last, 1]
    b <- axis[last, 2]
  }

  res$axes <- paste0("Parameters of the major axis (last iteration): a = ", a, ", b = ", b)
  res$puri <- paste("Item purification:", ifelse("purType" %in% names(model), model$purType, "unutilized"))

  res$thr <-
    paste0(
      "Detection threshold (last iteration): ",
      round(thr, 3),
      ifelse(thr == 1.5, " (fixed)", " (normal aproximation)")
    )

  res
})

output$DIF_DP_table_note <- renderUI({
  HTML(
    paste(
      "Notes:",
      DIF_DP_table_note()$puri,
      DIF_DP_table_note()$axes,
      DIF_DP_table_note()$thr,
      "Items detected as DIF are flagged with '***'",
      sep = "</br>"
    )
  )
})

# ** DIF items detected text ####
output$DIF_DP_dif_items <- renderPrint({
  DIFitems <- DIF_DP()$DIFitems
  if (DIFitems[1] == "no DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Purification info - number of iterations ####
output$DIF_DP_purification_info <- renderPrint({
  model <- DIF_DP()
  if (input$DIF_DP_purification & !is.null(DIF_DP()$difPur)) {
    cat("The table below describes the purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step, while the
        value of '0' means that the item was not detected as DIF. The first row corresponds to the initial classification of the items when all items
        were used for calculation of the DIF matching criterion. ")
    nrIter <- model$nrIter - 1
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_DP_purification & is.null(DIF_DP()$difPur)) {
    cat("No DIF items were detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_DP_purification_table <- reactive({
  model <- DIF_DP()
  tab <- model$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_DP_purification_table <- renderTable(
  {
    DIF_DP_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_DP_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_DP_purification.csv"
  },
  content = function(file) {
    data <- DIF_DP_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** REPORTS ####

# ** Delta plot calculations for report ####
report_DIF_DP <- reactive({
  data <- data.table(binary(), group = group())
  if (!input$customizeCheck) {
    DIF_DP_threshold_report <- input$DIF_DP_threshold
    purify_report <- input$DIF_DP_purification
    purType_report <- input$DIF_DP_purification_type
  } else {
    DIF_DP_threshold_report <- input$DIF_DP_threshold_report
    purify_report <- input$DIF_DP_purification_report
    purType_report <- input$DIF_DP_purification_type_report
  }

  switch(DIF_DP_threshold_report,
    "Fixed" = deltaPlot(data,
      group = "group",
      focal.name = 1,
      thr = 1.5,
      purify = purify_report
      # purType = purType_report
    ),
    "Normal" = deltaPlot(data,
      group = "group",
      focal.name = 1,
      thr = "norm",
      purify = purify_report,
      purType = purType_report
    )
  )
})

# ** Delta plot for report ####
report_DIF_DP_plot <- reactive({
  dp <- report_DIF_DP()
  df <- data.frame(dp$Deltas)
  df$nam <- item_numbers()

  par <- dp$axis.par
  thr <- dp$thr

  if (length(par) > 2) {
    par <- par[length(par) / 2, ]
  }

  if (length(thr) > 1) {
    thr <- thr[length(thr)]
  }

  p <- ggplot(
    df,
    aes(x = X1, y = X2, label = nam)
  ) +
    geom_point() +
    geom_text(hjust = 0, nudge_x = 0.05) +
    geom_abline(
      intercept = par[1], slope = par[2],
      linewidth = 1
    ) +
    geom_abline(
      intercept = par[1] + thr * sqrt(par[2]^2 + 1),
      slope = par[2],
      color = "red",
      linetype = "dashed",
      linewidth = 1
    ) +
    geom_abline(
      intercept = par[1] - thr * sqrt(par[2]^2 + 1),
      slope = par[2],
      color = "red",
      linetype = "dashed",
      linewidth = 1
    ) +
    labs(
      x = "Reference group",
      y = "Focal group"
    ) +
    xlim(min(dp$Deltas, na.rm = TRUE) - 0.5, max(dp$Deltas, na.rm = TRUE) + 0.5) +
    ylim(min(dp$Deltas, na.rm = TRUE) - 0.5, max(dp$Deltas, na.rm = TRUE) + 0.5) +
    theme_app()

  if (is.numeric(dp$DIFitems)) {
    df2 <- df[dp$DIFitems, ]
    p <- p + geom_point(
      data = df2,
      aes(x = X1, y = X2, label = nam),
      size = 6, color = "black", shape = 1
    )
  }
  p <- p + ggtitle("Delta plot")
  p
})


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * MANTEL-HAENSZEL ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

# ** Enabling/disabling options for combination of purification and correction ####
DIF_MH_combination_enable <- reactive(
  (input$DIF_MH_summary_purification & input$DIF_MH_summary_correction != "none")
)

observeEvent(DIF_MH_combination_enable(), {

  options <- if (DIF_MH_combination_enable()) {
    c("simple" = "simple", "iterative" = "iterative")
  } else {
    c("not applied" = "simple")
  }

  # updating selected choices for combinations
  updateSelectInput(
    session = session,
    label = "Combination",
    inputId = "DIF_MH_summary_combination",
    choices = options,
    selected = "simple"
  )
})

# ** Updating item and score sliders ####
observe({
  item_count <- ncol(binary())

  updateSliderInput(
    session = session,
    inputId = "DIF_MH_items_item",
    max = item_count
  )
  updateSliderInput(
    session = session,
    inputId = "DIF_MH_items_score",
    max = item_count,
    value = median(total_score(), na.rm = TRUE)
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MODEL ####

# ** Model ####
DIF_MH_model <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())

  fit <- ShinyItemAnalysis:::.difMH_edited(
    Data = data, group = group, focal.name = 1,
    p.adjust.method = input$DIF_MH_summary_correction,
    purify = input$DIF_MH_summary_purification, puriadjType = input$DIF_MH_summary_combination
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY ####

# ** DIF items detected text ####
output$DIF_MH_dif_items <- renderPrint({
  DIFitems <- DIF_MH_model()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_MH_summary_table <- reactive({
  res <- DIF_MH_model()
  deltaMH <- -2.35 * log(res$alphaMH)
  effsize <- symnum(abs(deltaMH), c(0, 1, 1.5, Inf), symbols = c("A", "B", "C"))
  blank <- character(length(res$alphaMH))

  pval <- if (res$p.adjust.method == "none") {
    round(res$p.value, 2)
  } else {
    round(res$adjusted.p, 2)
  }
  symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )

  tab <-
    data.frame(
      round(res$MH, 2),
      pval,
      symb,
      blank,
      res$alphaMH,
      deltaMH,
      effsize
    )

  colnames(tab) <-
    c(
      "MH (\\(\\mathit{\\chi^2}\\))",
      ifelse(
        res$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      "",
      "",
      "\\(\\mathit{\\alpha}_{\\mathrm{MH}}\\)",
      "\\(\\mathit{\\Delta}_{\\mathrm{MH}}\\)",
      ""
    )

  rownames(tab) <- item_names()

  tab
})

output$DIF_MH_summary_coef <- renderTable(
  {
    DIF_MH_summary_table()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Download summary table ####
output$DIF_MH_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_MH_statistics.csv"
  },
  content = function(file) {
    data <- DIF_MH_summary_table()
    data <- data[, -4]

    if ("adj. \\(\\mathit{p}\\)-value" %in% colnames(data)) {
      colnames(data) <- c("MH (X^2)", "adj. p-value", "sig. symb.", "MH alpha", "MH delta", "eff. size")
    } else {
      colnames(data) <- c("MH (X^2)", "p-value", "sig. symb.", "MH alpha", "MH delta", "eff. size")
    }

    write.csv(data, file)
    write(paste(
      "Note:",
      DIF_MH_summary_table_note()$p_adj,
      DIF_MH_summary_table_note()$puri,
      DIF_MH_summary_table_note()$combination,
      DIF_MH_summary_table_note()$thr_rounded,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      "DIF effect sizes (ETS Delta scale):
      A = negligible effect; B = moderate effect; C = large effect",
      "Effect size codes (for absolute values of 'MH delta'): 0 'A' 1.0 'B' 1.5 'C'",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Note about setting below summary table ####
DIF_MH_summary_table_note <- reactive({
  res <- NULL

  model <- DIF_MH_model()
  thr <- DIF_MH_model()$thr

  res$p_adj <- paste("P-value correction method:", switch(model$p.adjust.method,
    "BH" = "Benjamini-Hochberg",
    "BY" = "Benjamini-Yekutieli",
    "bonferroni" = "Bonferroni",
    "holm" = "Holm",
    "hochberg" = "Hochberg",
    "hommel" = "Hommel",
    "none" = "none"
  ))
  res$puri <- paste("Item purification:", ifelse(model$purification, "used", "unutilized"))
  res$combination <- ifelse(model$purification & model$p.adjust.method != "none",
                            paste("Item purification and p-value correction used in", input$DIF_MH_summary_combination, "combination"),
                            "")
  res$thr_rounded <- paste("Detection threshold:", round(thr[length(thr)], 3))

  res
})

output$DIF_MH_summary_table_note <- renderUI({
  HTML(
    paste(
      "Notes:",
      DIF_MH_summary_table_note()$p_adj,
      DIF_MH_summary_table_note()$puri,
      DIF_MH_summary_table_note()$combination,
      DIF_MH_summary_table_note()$thr_rounded,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      "DIF effect sizes (ETS Delta scale):
      <i>A</i> = negligible effect; <i>B</i> = moderate effect; <i>C</i> = large effect<br>
      Effect size codes (for absolute values of 'MH delta'): 0 'A' 1.0 'B' 1.5 'C'",
      sep = "</br>"
    )
  )
})

# ** Purification info - number of iterations ####
output$DIF_MH_summary_purification_info <- renderPrint({
  model <- DIF_MH_model()
  if (input$DIF_MH_summary_purification & !is.null(DIF_MH_model()$difPur)) {
    cat("The table below describes the purification process. The rows correspond to the purification iteration and the
        columns correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_MH_summary_purification & is.null(DIF_MH_model()$difPur)) {
    cat("No DIF items were detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_MH_summary_purification_table <- reactive({
  model <- DIF_MH_model()
  tab <- model$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_MH_summary_purification_table <- renderTable(
  {
    DIF_MH_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_MH_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_MH_purification.csv"
  },
  content = function(file) {
    data <- DIF_MH_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Contingency tables ####
DIF_MH_items_table <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())
  total <- total_score()
  item <- input$DIF_MH_items_item
  score <- input$DIF_MH_items_score

  df <- data.frame("Answer" = data[, item], "Group" = group)
  df$Answer <- relevel(
    factor(df$Answer, labels = c("Incorrect", "Correct")),
    "Correct"
  )
  df$Group <- factor(df$Group, labels = c("Reference group", "Focal group"))

  df <- df[total == score, ]

  tab <- xtabs(~ Group + Answer, data = df) |>
    addmargins(FUN = list(Total = sum), quiet = TRUE) |>
    as.data.frame.matrix()

  tab
})

# ** Contingency tables output ####
output$DIF_MH_items_table <- renderTable(
  {
    DIF_MH_items_table()
  },
  rownames = TRUE,
  digits = 0
)

# ** OR calculation ####
output$DIF_MH_items_interpretation <- renderUI({
  tab <- DIF_MH_items_table()
  a <- tab[1, 1]
  b <- tab[1, 2]
  c <- tab[2, 1]
  d <- tab[2, 2]
  OR <- round((a * d) / (b * c), 2)

  item <- input$DIF_MH_items_item
  score <- input$DIF_MH_items_score

  alphaMH <- round(DIF_MH_model()$alphaMH[item], 2)
  deltaMH <- -2.35 * log(alphaMH)
  effect_size <- symnum(abs(deltaMH), cutpoints = c(0, 1, 1.5, Inf), symbols = LETTERS[1:3])

  validate(
    need(((b * c > 0) & (a * d > 0)), "Odds ratio cannot be calculated!"),
    need((a > 0 & b > 0), "Not enough complete observations in the reference group."),
    need((c > 0 & d > 0), "Not enough complete observations in the focal group."),
    errorClass = "validation-error"
  )

  txt <- paste0(
    "For a respondent who reached the total score of ", score,
    ", the odds of answering item ", item_numbers()[item],
    " correctly is ",
    ifelse(OR == 1, " the same for both groups. ",
      ifelse(OR > 1,
        paste0(OR, " times higher in the reference group than in the focal group."),
        paste0(OR, " times lower in the reference group than in the focal group.")
      )
    )
  )

  txtMH <- paste0(
    "The Mantel-Haenszel estimate of the odds ratio accounting for all levels of the total score is equal to ",
    alphaMH, ". The odds of answering item ", item_numbers()[item],
    " correctly is ",
    ifelse(alphaMH == 1, " the same for both groups. ",
      ifelse(alphaMH > 1,
        paste0(alphaMH, " times higher in the reference group than in the focal group."),
        paste0(alphaMH, " times lower in the reference group than in the focal group.")
      )
    )
  )
  txtDelta <- paste0(
    "Mantel-Haenszel D-DIF index is equal to ", round(deltaMH, 2), ". This indicates category ", effect_size,
    " DIF effect size - ",
    switch(effect_size,
      "A" = "negligible",
      "B" = "moderate",
      "C" = "large"
    ), "."
  )
  HTML(paste(
    sprintf(
      paste("$$\\mathrm{OR} = \\frac{%d \\cdot %d}{%d \\cdot %d} = %.2f$$", txt),
      a, d, b, c, OR
    ), "<br/><br/>",
    sprintf(
      paste("$$\\alpha_{\\mathrm{MH}} = %.2f$$", txtMH),
      alphaMH
    ), "<br/><br/>",
    sprintf(
      paste("$$\\Delta_{\\mathrm{MH}} = -2.35 \\cdot \\log(\\alpha_{\\mathrm{MH}}) = %.2f$$", txtDelta),
      deltaMH
    )
  ))
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** REPORT ####

# ** Model for report ####
report_DIF_MH_model <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())

  if (!input$customizeCheck) {
    fit <- DIF_MH_model()
  } else {
    p.adjust.method_report <- input$report_DIF_MH_correction_method
    purify_report <- input$report_DIF_MH_purification

    fit <- ShinyItemAnalysis:::.difMH_edited(
      Data = data, group = group, focal.name = 1,
      p.adjust.method = p.adjust.method_report,
      purify = purify_report
    )
  }
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * SIBTEST ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

# ** Enabling/disabling options for combination of purification and correction ####
DIF_SIBTEST_combination_enable <- reactive(
  (input$DIF_SIBTEST_purification & input$DIF_SIBTEST_correction != "none")
)

observeEvent(DIF_SIBTEST_combination_enable(), {

  options <- if (DIF_SIBTEST_combination_enable()) {
    c("simple" = "simple", "iterative" = "iterative")
  } else {
    c("not applied" = "simple")
  }

  # updating selected choices for combinations
  updateSelectInput(
    session = session,
    label = "Combination",
    inputId = "DIF_SIBTEST_combination",
    choices = options,
    selected = "simple"
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_SIBTEST_model <- reactive({
  # data
  group <- unlist(group())
  data <- data.frame(binary())

  data_group <- bind_cols(data, group = group)
  group_counts <- data_group |>
    drop_na() |>
    count(group)
  sc0 <- group_counts |>
    filter(group == 0) |>
    pull(n)
  sc1 <- group_counts |>
    filter(group == 1) |>
    pull(n)

  validate(
    need(sc0 > 1L, "Not enough complete observations in the reference group."),
    need(sc1 > 1L, "Not enough complete observations in the focal group."),
    errorClass = "validation-error"
  )

  fit <- ShinyItemAnalysis:::.difSIBTEST_edited(
    Data = data, group = group, focal.name = 1,
    type = input$DIF_SIBTEST_type,
    purify = input$DIF_SIBTEST_purification,
    p.adjust.method = input$DIF_SIBTEST_correction,
    puriadjType = input$DIF_SIBTEST_combination
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MAIN TAB ####

# ** Warning for missing values ####
output$DIF_SIBTEST_NA_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_SIBTEST_dif_items <- renderPrint({
  DIFitems <- DIF_SIBTEST_model()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_SIBTEST_table <- reactive({
  res <- DIF_SIBTEST_model()

  # deal with only one pval base od model specs
  pval <- if (res$p.adjust.method == "none") {
    res$p.value
  } else {
    res$adjusted.p
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )

  tab <- data.frame(
    res$Beta,
    res$SE,
    res$X2,
    pval,
    pval_symb
  )

  colnames(tab) <-
    c(
      "\\(\\mathit{\\beta}\\)",
      "SE(\\(\\mathit{\\beta}\\))",
      "\\(\\mathit{\\chi^2}\\)",
      ifelse(
        res$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      ""
    )

  rownames(tab) <- item_names()

  tab
})

output$DIF_SIBTEST_coef <- renderTable(
  {
    DIF_SIBTEST_table()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_SIBTEST_table_note <- reactive({
  res <- NULL

  model <- DIF_SIBTEST_model()

  res$type <- paste0("Tested DIF type: ", switch(model$type,
    "udif" = "uniform",
    "nudif" = "non-uniform"
  ))

  res$p_adj <-
    paste("P-value correction method:", switch(model$p.adjust.method,
      bonferroni = "Bonferroni",
      holm = "Holm",
      hochberg = "Hochberg",
      hommel = "Hommel",
      BH = "Benjamini-Hochberg",
      BY = "Benjamini-Yekutieli",
      none = "none"
    ))

  res$puri <- paste("Item purification:", ifelse(model$purification, "used", "unutilized"))
  res$combination <- ifelse(model$purification & model$p.adjust.method != "none",
                            paste("Item purification and p-value correction used in", input$DIF_SIBTEST_combination, "combination"),
                            "")

  thr <- qchisq(1 - model$alpha, model$df[1])
  res$thr <-
    paste0(
      "Detection threshold: ",
      round(thr, 4),
      " (significance level: ",
      model$alpha,
      ")"
    )

  res
})

output$DIF_SIBTEST_table_note <- renderUI({
  HTML(
    paste(
      "Notes:",
      DIF_SIBTEST_table_note()$type,
      DIF_SIBTEST_table_note()$p_adj,
      DIF_SIBTEST_table_note()$puri,
      DIF_SIBTEST_table_note()$combination,
      DIF_SIBTEST_table_note()$thr,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_SIBTEST_table_download <- downloadHandler(
  filename = function() {
    paste("DIF_SIBTEST_statistics", ".csv", sep = "")
  },
  content = function(file) {
    data <- DIF_SIBTEST_table()

    colnames(data) <-
      c(
        "beta",
        "SE(beta)",
        "X^2",
        ifelse(
          "\\(\\mathit{p}\\)-value" %in% colnames(data),
          "p-value",
          "adj. p-value"
        ),
        "sig. symb."
      )

    rownames(data) <- item_names()

    write.csv(data, file)
    write(paste(
      "Notes:",
      DIF_SIBTEST_table_note()$type,
      DIF_SIBTEST_table_note()$p_adj,
      DIF_SIBTEST_table_note()$puri,
      DIF_SIBTEST_table_note()$combination,
      DIF_SIBTEST_table_note()$thr,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_SIBTEST_purification_info <- renderPrint({
  model <- DIF_SIBTEST_model()
  if (input$DIF_SIBTEST_purification & !is.null(model$difPur)) {
    cat("The table below describes the purification process. The rows correspond to the purification iteration and the
        columns correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_SIBTEST_purification & is.null(model$difPur)) {
    cat("No DIF items were detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_SIBTEST_purification_table <- reactive({
  tab <- DIF_SIBTEST_model()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})
output$DIF_SIBTEST_purification_table <- renderTable(
  {
    DIF_SIBTEST_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_SIBTEST_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_SIBTEST_purification.csv"
  },
  content = function(file) {
    data <- DIF_SIBTEST_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * LOGISTIC ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

match_logistic <- c("DIF_logistic_summary_matching", "DIF_logistic_items_matching")
puri_logistic <- c("DIF_logistic_summary_purification", "DIF_logistic_items_purification")

# ** Updating DIF matching & disable purification if DMV present ####
observe({
  if (DIF_matching_presence()) {
    lapply(match_logistic, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore",
          "Uploaded" = "uploaded",
          "Standardized uploaded" = "zuploaded"
        ),
        selected = "zscore"
      )
    })
  } else {
    lapply(match_logistic, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore"
        ),
        selected = "zscore"
      )
    })
  }
})

mapply(
  function(match, puri) {
    observeEvent(input[[paste0(match)]], {
      if (input[[paste0(match)]] %in% c("uploaded", "zuploaded")) {
        updateCheckboxInput(session, paste0(puri), value = FALSE)
        shinyjs::disable(paste0(puri))
      } else {
        shinyjs::enable(paste0(puri))
      }
    })
  },
  match = match_logistic, puri = puri_logistic
)

DIF_logistic <- reactiveValues(
  type = NULL,
  matching = NULL,
  parametrization = NULL,
  correction = NULL,
  purification = NULL,
  combination = NULL
)

# ** Updating type ####
observeEvent(input$DIF_logistic_summary_type, {
  DIF_logistic$type <- input$DIF_logistic_summary_type
})
observeEvent(input$DIF_logistic_items_type, {
  DIF_logistic$type <- input$DIF_logistic_items_type
})
observeEvent(DIF_logistic$type, {
  if (DIF_logistic$type != input$DIF_logistic_summary_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_logistic_summary_type",
      selected = DIF_logistic$type
    )
  }
  if (DIF_logistic$type != input$DIF_logistic_items_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_logistic_items_type",
      selected = DIF_logistic$type
    )
  }
})

# ** Updating matching ####
observeEvent(input$DIF_logistic_summary_matching, {
  DIF_logistic$matching <- input$DIF_logistic_summary_matching
})
observeEvent(input$DIF_logistic_items_matching, {
  DIF_logistic$matching <- input$DIF_logistic_items_matching
})
observeEvent(DIF_logistic$matching, {
  if (DIF_logistic$matching != input$DIF_logistic_summary_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_logistic_summary_matching",
      value = DIF_logistic$matching
    )
  }
  if (DIF_logistic$matching != input$DIF_logistic_items_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_logistic_items_matching",
      value = DIF_logistic$matching
    )
  }
})

# ** Updating parametrization ####
observeEvent(input$DIF_logistic_summary_parametrization, {
  DIF_logistic$parametrization <- input$DIF_logistic_summary_parametrization
})
observeEvent(input$DIF_logistic_items_parametrization, {
  DIF_logistic$parametrization <- input$DIF_logistic_items_parametrization
})
observeEvent(DIF_logistic$parametrization, {
  if (DIF_logistic$parametrization != input$DIF_logistic_summary_parametrization) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_logistic_summary_parametrization",
      selected = DIF_logistic$parametrization
    )
  }
  if (DIF_logistic$parametrization != input$DIF_logistic_items_parametrization) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_logistic_items_parametrization",
      selected = DIF_logistic$parametrization
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_logistic_summary_correction, {
  DIF_logistic$correction <- input$DIF_logistic_summary_correction
})
observeEvent(input$DIF_logistic_items_correction, {
  DIF_logistic$correction <- input$DIF_logistic_items_correction
})
observeEvent(DIF_logistic$correction, {
  if (DIF_logistic$correction != input$DIF_logistic_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_logistic_summary_correction",
      selected = DIF_logistic$correction
    )
  }
  if (DIF_logistic$correction != input$DIF_logistic_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_logistic_items_correction",
      selected = DIF_logistic$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_logistic_summary_purification, {
  DIF_logistic$purification <- input$DIF_logistic_summary_purification
})
observeEvent(input$DIF_logistic_items_purification, {
  DIF_logistic$purification <- input$DIF_logistic_items_purification
})
observeEvent(DIF_logistic$purification, {
  if (DIF_logistic$purification != input$DIF_logistic_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_logistic_summary_purification",
      value = DIF_logistic$purification
    )
  }
  if (DIF_logistic$purification != input$DIF_logistic_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_logistic_items_purification",
      value = DIF_logistic$purification
    )
  }
})

# ** Enabling/disabling options for combination of purification and correction ####
DIF_logistic_combination_enable <- reactive(
  (input$DIF_logistic_summary_purification & input$DIF_logistic_summary_correction != "none") |
    (input$DIF_logistic_items_purification & input$DIF_logistic_items_correction != "none")
)

observeEvent(DIF_logistic_combination_enable(), {

  options <- if (DIF_logistic_combination_enable()) {
    c("simple" = "simple", "iterative" = "iterative")
  } else {
    c("not applied" = "simple")
  }

  # updating selected choices for combinations
  updateSelectInput(
    session = session,
    label = "Combination",
    inputId = "DIF_logistic_summary_combination",
    choices = options,
    selected = "simple"
  )
  updateSelectInput(
    session = session,
    label = "Combination",
    inputId = "DIF_logistic_items_combination",
    choices = options,
    selected = "simple"
  )
})

# ** Updating combination ####
observeEvent(input$DIF_logistic_summary_combination, {
  DIF_logistic$combination <- input$DIF_logistic_summary_combination
})
observeEvent(input$DIF_logistic_items_combination, {
  DIF_logistic$combination <- input$DIF_logistic_items_combination
})
observeEvent(DIF_logistic$combination, {
  if (DIF_logistic$combination != input$DIF_logistic_summary_combination) {
    updateSelectInput(
      session = session,
      inputId = "DIF_logistic_summary_combination",
      selected = DIF_logistic$combination
    )
  }
  if (DIF_logistic$combination != input$DIF_logistic_items_combination) {
    updateSelectInput(
      session = session,
      inputId = "DIF_logistic_items_combination",
      selected = DIF_logistic$combination
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(binary())
  updateSliderInput(
    session = session,
    inputId = "DIF_logistic_items",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_logistic_model <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())

  if (input$DIF_logistic_summary_matching == "uploaded") {
    match <- unlist(DIFmatching())
  } else if (input$DIF_logistic_summary_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else if (input$DIF_logistic_summary_matching == "zscore") {
    match <- "zscore"
  } else if (input$DIF_logistic_summary_matching == "score") {
    match <- "score"
  }

  fit <- tryCatch(
    ShinyItemAnalysis:::.difLogistic_edited(
      Data = data, group = group, match = match, focal.name = 1,
      type = input$DIF_logistic_summary_type,
      p.adjust.method = input$DIF_logistic_summary_correction,
      purify = input$DIF_logistic_summary_purification,
      puriadjType = input$DIF_logistic_summary_combination,
      all.cov = TRUE
    ),
    error = function(e) e
  )

  validate(
    need(
      class(fit) == "Logistic",
      paste0("This method cannot be used on this data. Error returned: ", fit$message)
    ),
    errorClass = "validation-error"
  )

  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY ####

# ** Matching variable in text ####
DIF_logistic_summary_matching_text <- reactive({
  txt1 <- if (input$DIF_logistic_summary_matching == "uploaded") {
    "X_p"
  } else if (input$DIF_logistic_summary_matching == "zuploaded") {
    "Z_p"
  } else if (input$DIF_logistic_summary_matching == "zscore") {
    "Z_p"
  } else if (input$DIF_logistic_summary_matching == "score") {
    "X_p"
  }
  txt1
})

output$DIF_logistic_summary_matching_text <- renderUI({
  HTML(paste0("<b>\\(", DIF_logistic_summary_matching_text(), "\\)</b>"))
})

# ** Equation ####
DIF_logistic_summary_equation <- reactive({
  txt1 <- DIF_logistic_summary_matching_text()

  txt2 <- if (input$DIF_logistic_summary_parametrization == "irt") {
    paste0(
      "(a_{i} + a_{i\\text{DIF}} G_p)",
      "(", txt1, " - b_{i} - b_{i\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0} + \\beta_{i1} ", txt1, " + \\beta_{i2} G_p + \\beta_{i3} ", txt1, ":G_p")
  }
  txt3 <- paste0(txt1, ", G_p")

  txt <- paste0("$$\\mathrm{P}(Y_{pi} = 1|", txt3, ") = \\frac{e^{", txt2, "}}{1 + e^{", txt2, "}}$$")
  txt
})

output$DIF_logistic_summary_equation <- renderUI({
  HTML(DIF_logistic_summary_equation())
})

# ** Warning for missing values ####
output$DIF_logistic_summary_NA_warning <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_logistic_summary_dif_items <- renderPrint({
  DIFitems <- DIF_logistic_model()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_logistic_summary_coef <- reactive({
  fit <- DIF_logistic_model()

  stat <- fit$Logistik
  pval <- if (fit$p.adjust.method == "none") {
    fit$p.value
  } else {
    fit$adjusted.p
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )
  blank <- character(length(fit$names))

  r2 <- fit$deltaR2
  zt <- symnum(r2, c(0, 0.13, 0.26, 1),
    symbols = c("A", "B", "C")
  )
  jg <- symnum(r2, c(0, 0.035, 0.07, 1),
    symbols = c("A", "B", "C")
  )

  tab_coef <- fit$logitPar
  list_vcov <- fit$cov.M1
  if (all(fit$DIFitems != "No DIF item detected")) {
    list_vcov[fit$DIFitems] <- fit$cov.M0[fit$DIFitems]
  }

  if (input$DIF_logistic_summary_parametrization == "irt") {
    tab_coef_irt <- data.frame(
      a = tab_coef[, 2],
      b = -tab_coef[, 1] / tab_coef[, 2],
      aDIF = tab_coef[, 4],
      bDIF = (tab_coef[, 1] * tab_coef[, 4] - tab_coef[, 2] * tab_coef[, 3]) / (tab_coef[, 2] * (tab_coef[, 2] + tab_coef[, 4]))
    )
    # SE using delta method
    list_vcov_irt <- lapply(1:nrow(tab_coef), function(i) {
      nams <- names(which(tab_coef[i, ] != 0))
      vcov_tmp <- matrix(
        0,
        nrow = 4, ncol = 4,
        dimnames = list(
          colnames = names(tab_coef[i, ]),
          rownames = names(tab_coef[i, ])
        )
      )
      vcov_tmp[nams, nams] <- list_vcov[[i]]

      ShinyItemAnalysis:::delta_ses(
        list(~x2, ~ -x1 / x2, ~x4, ~ (x1 * x4 - x2 * x3) / (x2 * (x2 + x4))),
        mean = tab_coef[i, ],
        cov = vcov_tmp
      )
    })
    tab_se_irt <- do.call(rbind, list_vcov_irt)
    tab_par <-
      cbind(tab_coef_irt, tab_se_irt)[, order(c(seq(ncol(tab_coef_irt)), seq(ncol(tab_se_irt))))]
    tab_par <- tab_par[, c(1:2, 5:6, 3:4, 7:8)]
  } else {
    tab_se <- fit$logitSe
    tab_par <-
      cbind(tab_coef, tab_se)[, order(c(seq(ncol(tab_coef)), seq(ncol(tab_se))))]
  }

  tab <- data.frame(
    formatC(stat, format = "f", digits = 3),
    # round is ignored by renderTable, but formatC locks decimals in downloadTable on the other hand...
    formatC(pval, format = "f", digits = 3),
    pval_symb,
    blank,
    formatC(r2, format = "f", digits = 3),
    zt,
    jg,
    blank,
    tab_par
  )

  colnames(tab)[9:16] <- if (input$DIF_logistic_summary_parametrization == "irt") {
    c(
      "\\(\\mathit{a}\\)",
      "SE(\\(\\mathit{a}\\))",
      "\\(\\mathit{a_{\\text{DIF}}}\\)",
      "SE(\\(\\mathit{a_{\\text{DIF}}}\\))",
      "\\(\\mathit{b}\\)",
      "SE(\\(\\mathit{b}\\))",
      "\\(\\mathit{b_{\\text{DIF}}}\\)",
      "SE(\\(\\mathit{b_{\\text{DIF}}}\\))"
    )
  } else {
    c(
      "\\(\\mathit{\\beta_0}\\)",
      "SE(\\(\\mathit{\\beta_0}\\))",
      "\\(\\mathit{\\beta_1}\\)",
      "SE(\\(\\mathit{\\beta_1}\\))",
      "\\(\\mathit{\\beta_2}\\)",
      "SE(\\(\\mathit{\\beta_2}\\))",
      "\\(\\mathit{\\beta_3}\\)",
      "SE(\\(\\mathit{\\beta_3}\\))"
    )
  }

  colnames(tab)[1:8] <-
    c(
      "LR (\\(\\mathit{\\chi^2}\\))",
      ifelse(fit$p.adjust.method == "none", "\\(\\mathit{p}\\)-value", "adj. \\(\\mathit{p}\\)-value"),
      "",
      "",
      "\\(\\mathit{R^2}\\)",
      "ZT",
      "JG",
      ""
    )
  rownames(tab) <- item_names()
  tab
})

output$DIF_logistic_summary_coef <- renderTable(
  {
    DIF_logistic_summary_coef()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_logistic_summary_table_note <- reactive({
  res <- NULL

  fit <- DIF_logistic_model()
  thr <- DIF_logistic_model()$thr

  res$matching <- paste("Observed score:", switch(fit$match,
    "score" = "total score",
    "zscore" = "standardized total score",
    "matching variable" = "uploaded"
  ))
  res$type <- paste("DIF type tested:", switch(fit$type,
    "both" = "any DIF ",
    "udif" = "uniform DIF ",
    "nudif" = "non-uniform DIF "
  ))
  res$correction <- paste("P-value correction method:", switch(fit$p.adjust.method,
    "BH" = "Benjamini-Hochberg",
    "BY" = "Benjamini-Yekutieli",
    "bonferroni" = "Bonferroni",
    "holm" = "Holm",
    "hochberg" = "Hochberg",
    "hommel" = "Hommel",
    "none" = "none"
  ))
  res$purification <- paste("Item purification:", ifelse(fit$purification, "used", "unutilized"))
  res$combination <- ifelse(fit$purification & fit$p.adjust.method != "none",
                            paste("Item purification and p-value correction used in",
                                  input$DIF_logistic_summary_combination, "combination"),
                            "")
  res$threshold <- paste("Detection threshold:", round(thr[length(thr)], 3))

  res
})

output$DIF_logistic_summary_table_note <- renderUI({
  note <- DIF_logistic_summary_table_note()

  HTML(
    paste(
      "Notes:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      note$combination,
      note$threshold,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      "DIF effect sizes (Nagelkerke's \\(R^2\\)): <i>A</i> = negligible effect, <i>B</i> = moderate effect, <i>C</i> = large effect",
      "Thresholds by Zumbo & Thomas (ZT): 0 'A' 0.13 'B' 0.26 'C' 1",
      "Thresholds by Jodoin & Gierl (JG): 0 'A' 0.035 'B' 0.07 'C' 1",
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_logistic_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_logistic_table.csv"
  },
  content = function(file) {
    table <- DIF_logistic_summary_coef()

    # some dirty patch
    colnames(table) <- colnames(table) |>
      str_remove_all("\\\\\\(\\\\mathit\\{") |>
      str_remove_all("\\}\\\\\\)") |>
      str_replace("\\{\\\\text\\{DIF\\}\\}", "DIF") |>
      str_replace("\\\\chi", "X") |>
      str_remove_all("[\\\\]")

    names(table)[3] <- "sig. symb."

    # remove blank cols
    table[, c(4, 8)] <- NULL

    note <- DIF_logistic_summary_table_note()

    write.csv(table, file)
    write(paste(
      "Note:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      note$combination,
      note$threshold,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      "DIF effect sizes (Nagelkerke's R^2): A = negligible effect, B = moderate effect, C = large effect",
      "Thresholds by Zumbo & Thomas (ZT): 0 'A' 0.13 'B' 0.26 'C' 1",
      "Thresholds by Jodoin & Gierl (JG): 0 'A' 0.035 'B' 0.07 'C' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_logistic_summary_purification_info <- renderPrint({
  model <- DIF_logistic_model()
  if (input$DIF_logistic_summary_purification & !is.null(DIF_logistic_model()$difPur)) {
    cat("The table below describes purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_logistic_summary_purification & is.null(DIF_logistic_model()$difPur)) {
    cat("No DIF item was detected whatsoever, nothing to show.")
  } else {
    cat("Item purification not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_logistic_summary_purification_table <- reactive({
  tab <- DIF_logistic_model()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_logistic_summary_purification_table <- renderTable(
  {
    DIF_logistic_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_logistic_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_logistic_purification.csv"
  },
  content = function(file) {
    table <- DIF_logistic_summary_purification_table()
    write.csv(table, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_logistic_items_NA_warning <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot ####
DIF_logistic_items_plot <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())
  item <- input$DIF_logistic_items

  if (input$DIF_logistic_items_matching == "uploaded") {
    match <- unlist(DIFmatching())
  } else if (input$DIF_logistic_items_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else if (input$DIF_logistic_items_matching == "zscore") {
    match <- "zscore"
  } else if (input$DIF_logistic_items_matching == "score") {
    match <- "score"
  }

  fit <- DIF_logistic_model()

  g <- plotDIFLogistic(fit,
    item = item, match = match, item.name = item_names()[item],
    Data = data, group = group
  )
  g
})

output$DIF_logistic_items_plot <- renderPlotly({
  g <- DIF_logistic_items_plot()
  p <- ggplotly(g)

  p$x$data[[1]]$text <- paste0(
    "Group: Reference", "<br />",
    "Match: ", p$x$data[[1]]$x, "<br />",
    "Probability: ", p$x$data[[1]]$y
  )
  p$x$data[[2]]$text <- paste0(
    "Group: Focal", "<br />",
    "Match: ", p$x$data[[2]]$x, "<br />",
    "Probability: ", p$x$data[[2]]$y
  )
  p$x$data[[3]]$text <- gsub("gr1", "Reference", p$x$data[[3]]$text)
  p$x$data[[3]]$text <- gsub("Probability", "Empirical probability", p$x$data[[3]]$text)
  p$x$data[[4]]$text <- gsub("gr2", "Focal", p$x$data[[4]]$text)
  p$x$data[[4]]$text <- gsub("Probability", "Empirical probability", p$x$data[[4]]$text)

  if (input$DIF_logistic_items_matching == "uploaded") {
    match <- "Uploaded match"
  } else if (input$DIF_logistic_items_matching == "zuploaded") {
    match <- "Uploaded standardized match"
  } else if (input$DIF_logistic_items_matching == "zscore") {
    match <- "Z-score"
  } else if (input$DIF_logistic_items_matching == "score") {
    match <- "Score"
  }

  for (i in 1:length(p$x$data)) {
    text <- gsub("Match", match, p$x$data[[i]]$text)
    text <- lapply(strsplit(text, split = "<br />"), unique)
    text <- unlist(lapply(text, paste, collapse = "<br />"))
    p$x$data[[i]]$text <- text
  }

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot ####
output$DIF_logistic_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_logistic_", item_names()[input$DIF_logistic_items], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_logistic_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Matching variable in text ####
output$DIF_logistic_items_matching_text <- renderUI({
  HTML(paste0("<b>\\(", DIF_logistic_summary_matching_text(), "\\)</b>"))
})

# ** Equation ####
output$DIF_logistic_items_equation <- renderUI({
  HTML(DIF_logistic_summary_equation())
})

# ** Table with coefficients ####
output$DIF_logistic_items_coef <- renderTable(
  {
    tab <- DIF_logistic_summary_coef()
    item <- input$DIF_logistic_items

    tab <- tab[item, -c(1:8)]
    tab_coef <- unlist(tab[, seq(1, ncol(tab), 2)])
    tab_se <- unlist(tab[, seq(2, ncol(tab), 2)])

    tab <- data.frame(tab_coef, tab_se)
    # rownames(tab) <- c("\\(\\mathit{b}_0\\)", "\\(\\mathit{b}_1\\)", "\\(\\mathit{b}_2\\)", "\\(\\mathit{b}_3\\)")
    colnames(tab) <- c("Estimate", "SE")

    tab
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** REPORTS ####

# ** Method fit for report ####
report_DIF_logistic_model <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())

  if (!input$customizeCheck) {
    fit <- DIF_logistic_model()
  } else {
    type_report <- input$type_print_DIF_logistic_report
    correction_report <- input$correction_method_log_report
    purify_report <- input$puri_LR_report

    fit <- ShinyItemAnalysis:::.difLogistic_edited(
      Data = data, group = group, focal.name = 1,
      type = type_report,
      p.adjust.method = correction_report,
      purify = purify_report
    )
  }
  fit
})

# ** Plot for report ####
report_DIF_logistic_plot <- reactive({
  group <- unlist(group())
  data <- data.frame(binary())

  fit <- report_DIF_logistic_model()

  if (fit$DIFitems[1] != "No DIF item detected") {
    graflist <- vector("list", length = length(fit$DIFitems))
    i <- 1
    for (item in fit$DIFitems) {
      g <- plotDIFLogistic(fit,
        item = item,
        item.name = item_names()[item],
        Data = data, group = group
      )
      g <- g + ggtitle(paste0("DIF logistic plot for ", item_names()[item])) +
        theme(
          text = element_text(size = 12),
          plot.title = element_text(size = 12, face = "bold")
        )
      graflist[[i]] <- g
      i <- i + 1
    }
  } else {
    graflist <- NULL
  }
  graflist
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * NLR ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

DIF_NLR <- reactiveValues(
  model = NULL,
  type = NULL,
  correction = NULL,
  purification = NULL,
  method = NULL
)

# ** Updating model ####
observeEvent(input$DIF_NLR_summary_model, {
  DIF_NLR$model <- input$DIF_NLR_summary_model
})
observeEvent(input$DIF_NLR_items_model, {
  DIF_NLR$model <- input$DIF_NLR_items_model
})
observeEvent(DIF_NLR$model, {
  if (DIF_NLR$model != input$DIF_NLR_summary_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_NLR_summary_model",
      selected = DIF_NLR$model
    )
  }
  if (DIF_NLR$model != input$DIF_NLR_items_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_NLR_items_model",
      selected = DIF_NLR$model
    )
  }
})

# ** Updating type ####
observeEvent(input$DIF_NLR_summary_type, {
  DIF_NLR$type <- input$DIF_NLR_summary_type
})
observeEvent(input$DIF_NLR_items_type, {
  DIF_NLR$type <- input$DIF_NLR_items_type
})
observeEvent(DIF_NLR$type, {
  if (length(DIF_NLR$type) != length(input$DIF_NLR_summary_type)) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_NLR_summary_type",
      selected = DIF_NLR$type
    )
  } else {
    if (any(DIF_NLR$type != input$DIF_NLR_summary_type)) {
      updateCheckboxGroupInput(
        session = session,
        inputId = "DIF_NLR_summary_type",
        selected = DIF_NLR$type
      )
    }
  }
  if (length(DIF_NLR$type) != length(input$DIF_NLR_items_type)) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_NLR_items_type",
      selected = DIF_NLR$type
    )
  } else {
    if (any(DIF_NLR$type != input$DIF_NLR_items_type)) {
      updateCheckboxGroupInput(
        session = session,
        inputId = "DIF_NLR_items_type",
        selected = DIF_NLR$type
      )
    }
  }
})

# ** Updating correction ####
observeEvent(input$DIF_NLR_summary_correction, {
  DIF_NLR$correction <- input$DIF_NLR_summary_correction
})
observeEvent(input$DIF_NLR_items_correction, {
  DIF_NLR$correction <- input$DIF_NLR_items_correction
})
observeEvent(DIF_NLR$correction, {
  if (DIF_NLR$correction != input$DIF_NLR_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_NLR_summary_correction",
      selected = DIF_NLR$correction
    )
  }
  if (DIF_NLR$correction != input$DIF_NLR_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_NLR_items_correction",
      selected = DIF_NLR$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_NLR_summary_purification, {
  DIF_NLR$purification <- input$DIF_NLR_summary_purification
})
observeEvent(input$DIF_NLR_items_purification, {
  DIF_NLR$purification <- input$DIF_NLR_items_purification
})
observeEvent(DIF_NLR$purification, {
  if (DIF_NLR$purification != input$DIF_NLR_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_summary_purification",
      value = DIF_NLR$purification
    )
  }
  if (DIF_NLR$purification != input$DIF_NLR_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_items_purification",
      value = DIF_NLR$purification
    )
  }
})

# ** Updating matching ####
observeEvent(input$DIF_NLR_summary_matching, {
  DIF_NLR$matching <- input$DIF_NLR_summary_matching
})
observeEvent(input$DIF_NLR_items_matching, {
  DIF_NLR$matching <- input$DIF_NLR_items_matching
})
observeEvent(DIF_NLR$matching, {
  if (DIF_NLR$matching != input$DIF_NLR_summary_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_summary_matching",
      value = DIF_NLR$matching
    )
  }
  if (DIF_NLR$matching != input$DIF_NLR_items_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_items_matching",
      value = DIF_NLR$matching
    )
  }
})

# ** Updating est. method ####
observeEvent(input$DIF_NLR_summary_method, {
  DIF_NLR$method <- input$DIF_NLR_summary_method
})
observeEvent(input$DIF_NLR_items_method, {
  DIF_NLR$method <- input$DIF_NLR_items_method
})
observeEvent(DIF_NLR$method, {
  if (DIF_NLR$method != input$DIF_NLR_summary_method) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_summary_method",
      value = DIF_NLR$method
    )
  }
  if (DIF_NLR$method != input$DIF_NLR_items_method) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_NLR_items_method",
      value = DIF_NLR$method
    )
  }
})

# create vectors with DMV and purification inputs names
# (to feed "lapplies" in respective method section)
match_NLR <- c("DIF_NLR_summary_matching", "DIF_NLR_items_matching")
puri_NLR <- c("DIF_NLR_summary_purification", "DIF_NLR_items_purification")

# update selectInput & disable purification if DMV present
observe({
  if (DIF_matching_presence()) {
    lapply(match_NLR, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Standardized total score" = "zscore",
          "Standardized uploaded" = "zuploaded"
        ),
        selected = "zscore"
      )
    })
  } else {
    lapply(match_NLR, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Standardized total score" = "zscore"
        ),
        selected = "zscore"
      )
    })
  }
})

mapply(
  function(match, puri) {
    observeEvent(input[[paste0(match)]], {
      if (input[[paste0(match)]] %in% c("uploaded", "zuploaded")) {
        updateCheckboxInput(session, paste0(puri), value = FALSE)
        shinyjs::disable(paste0(puri))
      } else {
        shinyjs::enable(paste0(puri))
      }
    })
  },
  match = match_NLR, puri = puri_NLR
)

# ** Updating item slider ####
observe({
  item_count <- ncol(binary())
  updateSliderInput(
    session = session,
    inputId = "DIF_NLR_items_item",
    max = item_count
  )
})

# ** Enabling/disabling options for type of DIF in summary tab ####
observeEvent(input$DIF_NLR_summary_model, {
  # what parameters can be selected with choice of model
  enaSelection <- switch(input$DIF_NLR_summary_model,
    "Rasch" = c("b"),
    "1PL" = c("b"),
    "2PL" = c("a", "b"),
    "3PLcg" = c("a", "b"),
    "3PLdg" = c("a", "b"),
    "3PLc" = c("a", "b", "c"),
    "3PLd" = c("a", "b", "d"),
    "4PLcgdg" = c("a", "b"),
    "4PLcg" = c("a", "b", "d"),
    "4PLdg" = c("a", "b", "c"),
    "4PL" = c("a", "b", "c", "d")
  )
  # what parameters cannot be selected with choice of model
  disSelection <- setdiff(letters[1:4], enaSelection)

  # converting letters to numbers
  myLetters <- letters[1:26]
  disNum <- match(disSelection, myLetters)
  enaNum <- match(enaSelection, myLetters)

  # updating selected choices for type of DIF
  updateCheckboxGroupInput(
    session = session,
    inputId = "DIF_NLR_summary_type",
    selected = enaSelection
  )

  # create object that identifies enabled and disabled options
  disElement <- paste0("#DIF_NLR_summary_type :nth-child(", disNum, ") label")
  enaElement <- paste0("#DIF_NLR_summary_type :nth-child(", enaNum, ") label")

  # disable checkbox options of group
  shinyjs::enable(selector = enaElement)
  shinyjs::disable(selector = disElement)
})

# ** Enabling/disabling options for type of DIF in items tab ####
observeEvent(input$DIF_NLR_items_model, {
  # what parameters can be selected with choice of model
  enaSelection <- switch(input$DIF_NLR_items_model,
    "Rasch" = c("b"),
    "1PL" = c("b"),
    "2PL" = c("a", "b"),
    "3PLcg" = c("a", "b"),
    "3PLdg" = c("a", "b"),
    "3PLc" = c("a", "b", "c"),
    "3PLd" = c("a", "b", "d"),
    "4PLcgdg" = c("a", "b"),
    "4PLcg" = c("a", "b", "d"),
    "4PLdg" = c("a", "b", "c"),
    "4PL" = c("a", "b", "c", "d")
  )
  # what parameters cannot be selected with choice of model
  disSelection <- setdiff(letters[1:4], enaSelection)

  # converting letters to numbers
  myLetters <- letters[1:26]
  disNum <- match(disSelection, myLetters)
  enaNum <- match(enaSelection, myLetters)

  # updating selected choices for type of DIF
  updateCheckboxGroupInput(
    session = session,
    inputId = "DIF_NLR_items_type",
    selected = enaSelection
  )

  # create object that identifies enabled and disabled options
  disElement <- paste0("#DIF_NLR_items_type :nth-child(", disNum, ") label")
  enaElement <- paste0("#DIF_NLR_items_type :nth-child(", enaNum, ") label")

  # disable checkbox options of group
  shinyjs::enable(selector = enaElement)
  shinyjs::disable(selector = disElement)
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_NLR_method <- reactive({
  data <- data.frame(binary())
  group <- unlist(group())

  model <- input$DIF_NLR_summary_model
  type <- paste0(input$DIF_NLR_summary_type, collapse = "")
  adj.method <- input$DIF_NLR_summary_correction
  purify <- input$DIF_NLR_summary_purification
  est_method <- input$DIF_NLR_summary_method

  # if (input$DIF_NLR_summary_matching == "zscore") {
  #   match <- z_score()
  # } else
  if (input$DIF_NLR_summary_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else {
    match <- input$DIF_NLR_summary_matching
  }

  fit <- tryCatch(
    difNLR(
      Data = data, group = group, focal.name = 1, match = match,
      model = model, type = type,
      p.adjust.method = adj.method, purify = purify,
      test = "LR", method = est_method
    ),
    error = function(e) e
  )

  validate(
    need(
      inherits(fit, "difNLR"),
      paste0("This method cannot be used on this data. Error returned: ", fit$message)
    ),
    errorClass = "validation-error"
  )

  fit
}) |>
  bindEvent(
    input$DIF_NLR_summary_run,
    input$DIF_NLR_items_run,
    binary(),
    group(),
    ignoreNULL = FALSE # ignore BTN for the first time, but require press for any further eval.
    )


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY TAB ####

# ** Equation ####
output$DIF_NLR_summary_equation <- renderUI({
  model <- input$DIF_NLR_summary_model

  if (model == "Rasch") {
    txta <- ""
  } else {
    if (model == "1PL") {
      txta <- "a_i"
    } else {
      txta <- "a_{iG_p}"
    }
  }

  txtb <- "b_{iG_p}"

  txt2 <- paste0(txta, "\\left(Z_p - ", txtb, "\\right)")
  txt2 <- paste0("e^{", txt2, "}")
  txt2 <- paste0("\\frac{", txt2, "}{1 + ", txt2, "}")

  if (model %in% c("3PLcg", "4PLcgdg", "4PLcg")) {
    txtc <- "c_i"
  } else {
    if (model %in% c("3PLc", "4PLdg", "4PL")) {
      txtc <- "c_{iG_p}"
    } else {
      txtc <- ""
    }
  }

  if (model %in% c("3PLdg", "4PLcgdg", "4PLdg")) {
    txtd <- "d_i"
  } else {
    if (model %in% c("3PLd", "4PLcg", "4PL")) {
      txtd <- "d_{iG_p}"
    } else {
      txtd <- ""
    }
  }

  if (txtc == "" & txtd == "") {
    txt3 <- ""
  } else {
    if (txtd == "") {
      txt3 <- paste0(txtc, " + \\left(1 - ", txtc, "\\right) \\cdot ")
    } else {
      if (txtc == "") {
        txt3 <- txtd
      } else {
        txt3 <- paste0(txtc, " + \\left(", txtd, " - ", txtc, "\\right) \\cdot ")
      }
    }
  }

  txt1 <- paste0(
    "\\mathrm{P}\\left(Y_{pi} = 1 | Z_p, G_p\\right) = "
  )

  txt <- paste0("$$", txt1, txt3, txt2, "$$")
  txt
})

# ** Warning for missing values ####
output$DIF_NLR_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_NLR_summary_dif_items <- renderPrint({
  DIFitems <- DIF_NLR_method()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Coefs ####
DIF_NLR_summary_parameters <- reactive({
  model <- DIF_NLR_method()

  cfs <- coef(model, IRTpars = TRUE, SE = TRUE, CI = 0, simplify = TRUE)

  cfs <- cfs |>
    as_tibble(rownames = ".item") |>
    separate_wider_delim(.item, delim = " ", names = c(".item", ".type")) |>
    filter(.type %in% c("estimate", "SE")) # so CIs cannot creep into the table

  cfs |>
    pivot_wider(
      names_from = .type, values_from = -c(.item, .type),
      names_glue = "{.type}_{.value}"
    ) |>
    select(-.item)
})
# ** Summary table ####
DIF_NLR_summary_table <- reactive({
  model <- DIF_NLR_method()
  stat <- model$Sval

  # deal with only one pval base od model specs
  if (model$p.adjust.method == "none") {
    pval <- model$pval
    pval_name <- "pval"
  } else {
    pval <- model$adj.pval
    pval_name <- "pval_adj"
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )
  pval_symb[pval_symb == "?"] <- ""

  cfs <- DIF_NLR_summary_parameters()

  tibble(
    stat,
    "{pval_name}" := pval,
    pval_symb,
    blank_col = "",
    cfs
  )
})

output$DIF_NLR_summary_coef <- renderTable(
  {
    tbl <- DIF_NLR_summary_table()

    rnm_vec <- c(
      "LR (\\(\\mathit{\\chi^2}\\))" = "stat",
      "\\(\\mathit{p}\\)-value" = "pval",
      "adj. \\(\\mathit{p}\\)-value" = "pval_adj"
    )

    tbl <- tbl |>
      rename(any_of(rnm_vec)) |> # use any_of to pick only existing columns (picking the right pval)
      rename_with(~ str_replace(.x, "Dif", "_{DIF}")) |>
      rename_with(~ paste0(str_replace(.x, "estimate_", "\\\\(\\\\mathit{"), "}\\)"), .cols = starts_with("estimate_")) |>
      rename_with(~ paste0(str_replace(.x, "SE_", "SE(\\\\(\\\\mathit{"), "}\\))"), .cols = starts_with("SE_"))


    tbl <- tbl |>
      mutate(rowname = item_names()) |>
      column_to_rownames()

    # rm names, which is not legal in tibble
    colnames(tbl)[colnames(tbl) %in% c("pval_symb", "blank_col")] <- ""

    tbl
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_NLR_summary_table_note <- reactive({
  res <- NULL

  model <- DIF_NLR_method()

  df <- model$df
  if (length(dim(df)) == 2L) {
    df <- df[, 1L]
  }

  thr <- if (length(unique(df)) == 1) {
    unique(qchisq(1 - model$alpha, df))
  } else {
    NULL
  }

  res$mod <- paste("Model:", switch(unique(model$model),
    "Rasch" = "Rasch model",
    "1PL" = "1PL model",
    "2PL" = "2PL model",
    "3PL" = "3PL model",
    "3PLcg" = "3PL model with fixed guessing for groups",
    "3PLdg" = "3PL model with fixed inattention parameter for groups",
    "3PLc" = "3PL model",
    "3PLd" = "3PL model with inattention parameter",
    "4PLcgdg" = "4PL model with fixed guessing and inattention parameter for groups",
    "4PLcgd" = "4PL model with fixed guessing for groups",
    "4PLd" = "4PL model with fixed guessing for groups",
    "4PLcdg" = "4PL model with fixed inattention parameter for groups",
    "4PLc" = "4PL model with fixed inattention parameter for groups",
    "4PL" = "4PL model"
  ))

  res$dmv <- paste("Observed score:", switch(as.character(model$match[1]), # ensures number is recognize as unnamed element
    "score" = "total score",
    "zscore" = "standardized total score",
    "uploaded"
  ))

  res$type <-
    paste0(
      "DIF type tested: difference in parameters ",
      paste0(input$DIF_NLR_summary_type, collapse = ", ")
    )

  res$p_adj <- paste("P-value correction method:", switch(model$p.adjust.method,
    holm = "Holm",
    hochberg = "Hochberg",
    hommel = "Hommel",
    bonferroni = "Bonferroni",
    BH = "Benjamini-Hochberg",
    BY = "Benjamini-Yekutieli",
    fdr = "FDR",
    none = "none"
  ))

  res$puri <- paste("Item purification:", ifelse(model$purification == T, "used", "unutilized"))
  res$thr_rounded <- paste("Detection threshold:", round(thr, 3))

  res
})

output$DIF_NLR_summary_table_note <- renderUI({
  HTML(
    paste(
      "Notes:",
      DIF_NLR_summary_table_note()$dmv,
      DIF_NLR_summary_table_note()$mod,
      DIF_NLR_summary_table_note()$type,
      DIF_NLR_summary_table_note()$p_adj,
      DIF_NLR_summary_table_note()$puri,
      DIF_NLR_summary_table_note()$thr_rounded,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "</br>"
    )
  )
})

# ** Download table ####
output$DIF_NLR_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_NLR_statistics.csv"
  },
  content = function(file) {
    tbl <- DIF_NLR_summary_table()

    rnm_vec <- c(
      "LR (X^2)" = "stat",
      "p-value" = "pval",
      "adj. p-value" = "pval_adj",
      "sig. symb." = "pval_symb"
    )

    tbl <- tbl |>
      select(-blank_col) |>
      rename(any_of(rnm_vec)) |>
      rename_with(~ str_replace(.x, "Dif", "_DIF")) |>
      rename_with(~ str_remove(.x, "estimate_"), .cols = starts_with("estimate_")) |>
      rename_with(~ paste0(str_replace(.x, "SE_", "SE("), ")"), .cols = starts_with("SE_"))


    tbl <- tbl |>
      mutate(rowname = item_names()) |>
      column_to_rownames()

    write.csv(tbl, file) # w/o blank col
    write(paste(
      "\nNote:",
      DIF_NLR_summary_table_note()$dmv,
      DIF_NLR_summary_table_note()$mod,
      gsub(",", "", DIF_NLR_summary_table_note()$type), # get rid of the comma - it separates col in CSV
      DIF_NLR_summary_table_note()$p_adj,
      DIF_NLR_summary_table_note()$puri,
      DIF_NLR_summary_table_note()$thr_rounded,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_NLR_summary_purification_info <- renderPrint({
  model <- DIF_NLR_method()
  if (input$DIF_NLR_summary_purification & !is.null(DIF_NLR_method()$difPur)) {
    cat("The table below describes purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$conv.puri, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_NLR_summary_purification & is.null(DIF_NLR_method()$difPur)) {
    cat("No DIF item was detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_NLR_summary_purification_table <- reactive({
  tab <- DIF_NLR_method()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_NLR_summary_purification_table <- renderTable(
  {
    DIF_NLR_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_NLR_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_NLR_purification.csv"
  },
  content = function(file) {
    data <- DIF_NLR_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS TAB ####

# ** Warning for missing values ####
output$DIF_NLR_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot ####
DIF_NLR_items_plot <- reactive({
  fit <- DIF_NLR_method()
  item <- input$DIF_NLR_items_item

  g <- plot(fit, item = item)[[1]] +
    theme_app() +
    theme(
      legend.box.just = "top",
      legend.position = "inside",
      legend.position.inside = c(0.01, 0.98),
      legend.justification = c(0, 1),
      legend.key.width = unit(1, "cm"),
      legend.box = "horizontal"
    ) +
    ggtitle(item_names()[item])
  g
})

output$DIF_NLR_items_plot <- renderPlotly({
  g <- DIF_NLR_items_plot()
  p <- ggplotly(g)

  p$x$data[[1]]$text <- paste0(
    "Group: Reference", "<br />",
    "Match: ", p$x$data[[1]]$x, "<br />",
    "Probability: ", p$x$data[[1]]$y
  )
  p$x$data[[2]]$text <- paste0(
    "Group: Focal", "<br />",
    "Match: ", p$x$data[[2]]$x, "<br />",
    "Probability: ", p$x$data[[2]]$y
  )

  p$x$data[[3]]$text <- gsub("size", "Group: Reference<br />Count", p$x$data[[3]]$text)
  p$x$data[[3]]$text <- gsub("match", "Z-score", p$x$data[[3]]$text)
  p$x$data[[3]]$text <- gsub("prob", "Empirical probability", p$x$data[[3]]$text)
  p$x$data[[4]]$text <- gsub("size", "Group: Focal<br />Count", p$x$data[[4]]$text)
  p$x$data[[4]]$text <- gsub("match", "Z-score", p$x$data[[4]]$text)
  p$x$data[[4]]$text <- gsub("prob", "Empirical probability", p$x$data[[4]]$text)

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot ####
output$DIF_NLR_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIFNonlinear_", item_names()[input$DIF_NLR_items_plot], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_NLR_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Equation ####
output$DIF_NLR_items_equation <- renderUI({
  model <- input$DIF_NLR_items_model

  if (model == "Rasch") {
    txta <- ""
  } else {
    if (model == "1PL") {
      txta <- "a_i"
    } else {
      txta <- "a_{iG_p}"
    }
  }

  txtb <- "b_{iG_p}"

  txt2 <- paste0(txta, "\\left(Z_p - ", txtb, "\\right)")
  txt2 <- paste0("e^{", txt2, "}")
  txt2 <- paste0("\\frac{", txt2, "}{1 + ", txt2, "}")

  if (model %in% c("3PLcg", "4PLcgdg", "4PLcg")) {
    txtc <- "c_i"
  } else {
    if (model %in% c("3PLc", "4PLdg", "4PL")) {
      txtc <- "c_{iG_p}"
    } else {
      txtc <- ""
    }
  }

  if (model %in% c("3PLdg", "4PLcgdg", "4PLdg")) {
    txtd <- "d_i"
  } else {
    if (model %in% c("3PLd", "4PLcg", "4PL")) {
      txtd <- "d_{iG_p}"
    } else {
      txtd <- ""
    }
  }

  if (txtc == "" & txtd == "") {
    txt3 <- ""
  } else {
    if (txtd == "") {
      txt3 <- paste0(txtc, " + \\left(1 - ", txtc, "\\right) \\cdot ")
    } else {
      if (txtc == "") {
        txt3 <- txtd
      } else {
        txt3 <- paste0(txtc, " + \\left(", txtd, " - ", txtc, "\\right) \\cdot ")
      }
    }
  }

  txt1 <- paste0(
    "\\mathrm{P}\\left(Y_{pi} = 1 | Z_p, G_p\\right) = "
  )

  txt <- paste0("$$", txt1, txt3, txt2, "$$")
  txt
})


# ** Table of item parameters ####
output$DIF_NLR_items_coef <- renderTable(
  {
    cfs <- DIF_NLR_summary_parameters()
    item_cfs <- cfs[input$DIF_NLR_items_item, ]

    item_cfs <- item_cfs |>
      pivot_longer(everything(), names_sep = "_", names_to = c(".value", "par"))

    item_cfs <- item_cfs |> column_to_rownames("par")

    rownames(item_cfs) <- paste0("\\(\\mathit{", gsub("Dif", "_{Dif}", rownames(item_cfs)), "}\\)")
    colnames(item_cfs) <- c("Estimate", "SE")

    item_cfs
  },
  include.rownames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * IRT LORD ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

DIF_Lord <- reactiveValues(
  model = NULL,
  correction = NULL,
  purification = NULL
)

# ** Updating model ####
observeEvent(input$DIF_Lord_summary_model, {
  DIF_Lord$model <- input$DIF_Lord_summary_model
})
observeEvent(input$DIF_Lord_items_model, {
  DIF_Lord$model <- input$DIF_Lord_items_model
})
observeEvent(DIF_Lord$model, {
  if (DIF_Lord$model != input$DIF_Lord_summary_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_Lord_summary_model",
      selected = DIF_Lord$model
    )
  }
  if (DIF_Lord$model != input$DIF_Lord_items_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_Lord_items_model",
      selected = DIF_Lord$model
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_Lord_summary_correction, {
  DIF_Lord$correction <- input$DIF_Lord_summary_correction
})
observeEvent(input$DIF_Lord_items_correction, {
  DIF_Lord$correction <- input$DIF_Lord_items_correction
})
observeEvent(DIF_Lord$correction, {
  if (DIF_Lord$correction != input$DIF_Lord_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_Lord_summary_correction",
      selected = DIF_Lord$correction
    )
  }
  if (DIF_Lord$correction != input$DIF_Lord_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_Lord_items_correction",
      selected = DIF_Lord$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_Lord_summary_purification, {
  DIF_Lord$purification <- input$DIF_Lord_summary_purification
})
observeEvent(input$DIF_Lord_items_purification, {
  DIF_Lord$purification <- input$DIF_Lord_items_purification
})
observeEvent(DIF_Lord$purification, {
  if (DIF_Lord$purification != input$DIF_Lord_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_Lord_summary_purification",
      value = DIF_Lord$purification
    )
  }
  if (DIF_Lord$purification != input$DIF_Lord_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_Lord_items_purification",
      value = DIF_Lord$purification
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(binary())
  updateSliderInput(
    session = session,
    inputId = "DIF_Lord_items_item",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_Lord_method <- reactive({
  group <- unlist(group(), use.names = FALSE)
  data <- data.frame(binary())

  if (input$DIF_Lord_summary_model == "3PL") {
    guess <- itemPar3PL(data)[, 3]
  }

  fit <- tryCatch(
    switch(input$DIF_Lord_summary_model,
      "1PL" = ShinyItemAnalysis:::.difLord_edited(
        Data = data, group = group, focal.name = 1,
        model = "1PL",
        p.adjust.method = input$DIF_Lord_summary_correction,
        purify = input$DIF_Lord_summary_purification
      ),
      "2PL" = ShinyItemAnalysis:::.difLord_edited(
        Data = data, group = group, focal.name = 1,
        model = "2PL",
        p.adjust.method = input$DIF_Lord_summary_correction,
        purify = input$DIF_Lord_summary_purification
      ),
      "3PL" = ShinyItemAnalysis:::.difLord_edited(
        Data = data, group = group, focal.name = 1,
        model = "3PL", c = guess,
        p.adjust.method = input$DIF_Lord_summary_correction,
        purify = input$DIF_Lord_summary_purification
      )
    ),
    error = function(e) e
  )

  validate(
    need(
      class(fit) == "Lord",
      paste0("This method cannot be used on this data. Error returned: ", fit$message)
    ),
    errorClass = "validation-error"
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY TAB ####

# ** Interpretation ####
output$DIF_Lord_summary_interpretation <- renderUI({
  type <- input$DIF_Lord_items_model

  txt <- switch(type,
    "1PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(b_{iR}\\) and \\(b_{iF}\\)
                             are difficulties for the reference and the focal group for item \\(i\\). "),
    "2PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\). "),
    "3PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters  \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\).
                             Parameter \\(c_i\\) is a common guessing parameter for item \\(i\\). ")
  )
  HTML(txt)
})

# ** Equation ####
output$DIF_Lord_summary_equation <- renderUI({
  type <- input$DIF_Lord_items_model
  eqR <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{\\theta_p - b_{iR}}}
                              {1 + e^{\\theta_p - b_{iR} }}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iR}
                              \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$")
  )

  eqF <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{\\theta_p - b_{iF}}}
                              {1 + e^{\\theta_p - b_{iF}}}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iF}
                              \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$")
  )
  paste(eqR, eqF)
})

# ** Warning for missing values ####
output$DIF_Lord_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_Lord_summary_dif_items <- renderPrint({
  DIFitems <- DIF_Lord_method()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_Lord_summary_table <- reactive({
  res <- DIF_Lord_method()

  # deal with only one pval base od model specs
  pval <- if (res$p.adjust.method == "none") {
    res$p.value
  } else {
    res$adjusted.p
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )
  blank <- character(length(res$names))

  # IRT pars
  if (res$purification) {
    m <- nrow(res$itemParFinal) / 2

    mR <- res$itemParFinal[1:m, ]
    mF <- res$itemParFinal[(m + 1):(2 * m), ]
  } else {
    m <- nrow(res$itemParInit) / 2

    mR <- res$itemParInit[1:m, ]
    mF <- res$itemParInit[(m + 1):(2 * m), ]

    mF <- itemRescale(mR, mF) # rescaling!
  }

  par <- rbind(mR, mF)

  wh_coef <- switch(res$model,
    "1PL" = 1,
    "2PL" = 1:2,
    "3PL" = c(1, 2)
  )
  wh_sd <- switch(res$model,
    "1PL" = 2,
    "2PL" = 3:4,
    "3PL" = 3:4
  )

  R_indices <- seq(1, m)
  F_indices <- seq(m + 1, 2 * m)

  if (res$model == "3PL") {
    guess_both <- par[R_indices, 6]
  }

  coefs_R <- par[R_indices, wh_coef]
  coefs_F <- par[F_indices, wh_coef]

  if (res$model != "1PL") {
    colnames(coefs_R) <- paste0(colnames(coefs_R), "_R")
    colnames(coefs_F) <- paste0(colnames(coefs_F), "_F")
  }

  coefs <- cbind(coefs_R, coefs_F)

  if (res$model == "1PL") {
    colnames(coefs) <- c("b_R", "b_F")
  }

  se_R <- par[R_indices, wh_sd]
  se_F <- par[F_indices, wh_sd]

  if (res$model != "1PL") {
    colnames(se_R) <- paste0(colnames(se_R), "_R")
    colnames(se_F) <- paste0(colnames(se_F), "_F")
  }

  ses <- cbind(se_R, se_F)

  if (res$model == "1PL") {
    colnames(ses) <- c("se(b)_R", "se(b)_F")
  }

  pars_zigzag <-
    cbind(coefs, ses)[, order(c(seq(ncol(coefs)), seq(ncol(ses))))]
  if (res$model == "3PL") {
    pars_zigzag <- cbind(pars_zigzag, guess_both)
  }


  colnames(pars_zigzag) <- switch(res$model,
    "1PL" = c(
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))"
    ),
    "2PL" = c(
      "\\(\\mathit{a}_{R}\\)",
      "SE(\\(\\mathit{a}_{R}\\))",
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{a}_{F}\\)",
      "SE(\\(\\mathit{a}_{F}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))"
    ),
    "3PL" = c(
      "\\(\\mathit{a}_{R}\\)",
      "SE(\\(\\mathit{a}_{R}\\))",
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{a}_{F}\\)",
      "SE(\\(\\mathit{a}_{F}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))",
      "\\(\\mathit{c}\\)"
    )
  )

  # from print.Lord:
  if (is.null(res$anchor.names)) {
    itk <- 1:length(res$LordChi)
  } else {
    itk <- (1:length(res$LordChi))[!is.na(res$LordChi)]
  }

  if (res$model == "1PL") {
    if (res$purification & is.null(res$anchor.names)) {
      pars <- res$itemParFinal
    } else {
      pars <- res$itemParInit
    }
    J <- nrow(pars) / 2
    mR <- pars[1:J, 1]
    mF <- itemRescale(pars[1:J, ], pars[(J + 1):(2 * J), ])[, 1]
    rr1 <- round(mF - mR, 4)
    rr2 <- round(-2.35 * rr1, 4)
    symb1 <-
      symnum(abs(rr2), c(0, 1, 1.5, Inf), symbols = c("A", "B", "C"))
    matR2 <- cbind(rr1, rr2)[itk, ]
  }

  tab <- if (res$model == "1PL") {
    data.frame(
      res$LordChi,
      pval,
      pval_symb,
      blank,
      matR2[, 1],
      matR2[, 2],
      symb1, # eff size based on delta
      blank,
      pars_zigzag
    )
  } else {
    data.frame(
      res$LordChi,
      pval,
      pval_symb,
      blank,
      pars_zigzag
    )
  }
  colnames(tab) <- if (res$model != "1PL") {
    c(
      "Lord's \\(\\mathit{\\chi^2}\\)",
      ifelse(
        res$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      "",
      "",
      colnames(pars_zigzag)
    )
  } else {
    c(
      "Lord's \\(\\mathit{\\chi^2}\\)",
      ifelse(
        res$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      "",
      "",
      "mF-mR",
      "deltaLord",
      "",
      "",
      colnames(pars_zigzag)
    )
  }

  rownames(tab) <- item_names()

  tab
})

# ** Note about setting below summary table ####
DIF_Lord_summary_table_note <- reactive({
  res <- NULL
  model <- DIF_Lord_method()

  res$mod <- paste0("Model: ", model$model)
  res$p_adj <- paste("P-value correction method:", switch(model$p.adjust.method,
    holm = "Holm",
    hochberg = "Hochberg",
    hommel = "Hommel",
    bonferroni = "Bonferroni",
    BH = "Benjamini-Hochberg",
    BY = "Benjamini-Yekutieli",
    none = "none"
  ))
  res$puri <- paste("Item purification:", ifelse(model$purification, "used", "unutilized"))
  res$thr <- paste0("Detection threshold: ", round(model$thr, 3))

  res
})

output$DIF_Lord_summary_table_note <- renderUI({
  txt_effect <- ifelse(input$DIF_Lord_summary_model == "1PL",
    "Effect size codes: 'A': negligible effect; 'B': moderate effect; 'C': large effect </br>
                       Effect size codes: 0 'A' 1.0 'B' 1.5 'C' (for absolute values of 'deltaLord')",
    ""
  )
  HTML(
    paste(
      "Notes:",
      DIF_Lord_summary_table_note()$mod,
      DIF_Lord_summary_table_note()$p_adj,
      DIF_Lord_summary_table_note()$puri,
      DIF_Lord_summary_table_note()$thr,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      txt_effect,
      sep = "</br>"
    )
  )
})

output$DIF_Lord_summary_coef <- renderTable(
  {
    DIF_Lord_summary_table()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Download summary table ####
output$DIF_Lord_summary_table_download <- downloadHandler(
  filename = function() {
    paste("DIF_Lord_statistics", ".csv", sep = "")
  },
  content = function(file) {
    data <- DIF_Lord_summary_table()
    data <- if (input$DIF_Lord_summary_model == "1PL") {
      data[, -c(4, 8)]
    } else {
      data[, -4]
    }

    # some dirty patch
    colnames(data) <- colnames(data) |>
      str_remove_all("\\\\\\(\\\\mathit\\{") |>
      str_remove_all("\\}\\\\\\)") |>
      str_replace("\\}_\\{", "_") |>
      str_replace("\\\\chi", "X") |>
      str_remove_all("[\\\\]")

    names(data)[3] <- "sig. symb."

    if (input$DIF_Lord_summary_model == "1PL") {
      names(data)[6] <- "eff. size"
    }

    rownames(data) <- item_names()

    write.csv(data, file)


    txt_effect <- ifelse(input$DIF_Lord_summary_model == "1PL",
      c(
        "Effect size codes: 'A': negligible effect; 'B': moderate effect; 'C': large effect",
        "Effect size codes: 0 'A' 1.0 'B' 1.5 'C' (for absolute values of 'deltaLord')"
      ),
      ""
    )

    write(paste(
      "Notes:",
      DIF_Lord_summary_table_note()$mod,
      DIF_Lord_summary_table_note()$p_adj,
      DIF_Lord_summary_table_note()$puri,
      DIF_Lord_summary_table_note()$thr,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      txt_effect,
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_Lord_summary_purification_info <- renderPrint({
  model <- DIF_Lord_method()
  if (input$DIF_Lord_summary_purification & !is.null(DIF_Lord_method()$difPur)) {
    cat("The table below describes purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_Lord_summary_purification & is.null(DIF_Lord_method()$difPur)) {
    cat("No DIF item was detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_Lord_summary_purification_table <- reactive({
  tab <- DIF_Lord_method()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_Lord_summary_purification_table <- renderTable(
  {
    DIF_Lord_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

output$DIF_Lord_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_Lord_purification.csv"
  },
  content = function(file) {
    data <- DIF_Lord_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_Lord_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot ####
DIF_Lord_items_plot <- reactive({
  fitLord <- DIF_Lord_method()
  item <- input$DIF_Lord_items_item

  g <- plotDIFirt(
    parameters = fitLord$itemParInit,
    item = item,
    item.name = item_names()[item]
  )[[item]]
  g
})

output$DIF_Lord_items_plot <- renderPlotly({
  g <- DIF_Lord_items_plot()
  p <- ggplotly(g)

  p$x$data[[1]]$text <- paste0(
    "Group: Reference", "<br />",
    "Ability: ", p$x$data[[1]]$x, "<br />",
    "Probability: ", p$x$data[[1]]$y
  )
  p$x$data[[2]]$text <- paste0(
    "Group: Focal", "<br />",
    "Ability: ", p$x$data[[2]]$x, "<br />",
    "Probability: ", p$x$data[[2]]$y
  )

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot ####
output$DIF_Lord_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIFIRTLord_", item_names()[input$DIF_Lord_items_item], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_Lord_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Interpretation ####
output$DIF_Lord_items_interpretation <- renderUI({
  type <- input$DIF_Lord_items_model

  txt <- switch(type,
    "1PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(b_{iR}\\) and \\(b_{iF}\\)
                             are difficulties for the reference and the focal group for item \\(i\\). "),
    "2PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\). "),
    "3PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters  \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\).
                             Parameter \\(c_i\\) is a common guessing parameter for item \\(i\\). ")
  )
  HTML(txt)
})

# ** Equation ####
output$DIF_Lord_items_equation <- renderUI({
  type <- input$DIF_Lord_items_model
  eqR <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{\\theta_p - b_{iR}}}
                              {1 + e^{\\theta_p - b_{iR} }}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iR}
                              \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$")
  )

  eqF <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{\\theta_p - b_{iF}}}
                              {1 + e^{\\theta_p - b_{iF}}}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iF}
                              \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$")
  )
  paste(eqR, eqF)
})

# ** Table of item parameters ####
DIF_Lord_items_coef <- reactive({
  fitLord <- DIF_Lord_method()

  if (input$DIF_Lord_items_purification) {
    m <- nrow(fitLord$itemParFinal) / 2
    mR <- fitLord$itemParFinal[1:m, ]
    mF <- fitLord$itemParFinal[(m + 1):(2 * m), ]
  } else {
    m <- nrow(fitLord$itemParInit) / 2
    mR <- fitLord$itemParInit[1:m, ]
    mF <- fitLord$itemParInit[(m + 1):(2 * m), ]
    mF <- itemRescale(mR, mF)
  }

  par <- rbind(mR, mF)

  wh_coef <- switch(input$DIF_Lord_items_model,
    "1PL" = 1,
    "2PL" = 1:2,
    "3PL" = c(1, 2, 6)
  )
  wh_sd <- switch(input$DIF_Lord_items_model,
    "1PL" = 2,
    "2PL" = 3:4,
    "3PL" = 3:4
  )

  item <- input$DIF_Lord_items_item
  tab_coef <- c(par[c(item, m + item), wh_coef])
  tab_sd <- c(par[c(item, m + item), wh_sd])

  if (input$DIF_Lord_items_model == "3PL") {
    tab_coef <- tab_coef[-6]
  }

  if (input$DIF_Lord_items_model == "3PL") {
    tab_sd <- c(tab_sd, NA)
  }

  tab <- data.frame(tab_coef, tab_sd)
  rownames(tab) <- switch(input$DIF_Lord_items_model,
    "1PL" = c("\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)"),
    "2PL" = c("\\(\\mathit{a}_{R}\\)", "\\(\\mathit{a}_{F}\\)", "\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)"),
    "3PL" = c("\\(\\mathit{a}_{R}\\)", "\\(\\mathit{a}_{F}\\)", "\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)", "\\(\\mathit{c}\\)")
  )
  colnames(tab) <- c("Estimate", "SE")

  tab
})

output$DIF_Lord_items_coef <- renderTable(
  {
    DIF_Lord_items_coef()
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * IRT Raju ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

DIF_Raju <- reactiveValues(
  model = NULL,
  correction = NULL,
  purification = NULL
)

# ** Updating model ####
observeEvent(input$DIF_Raju_summary_model, {
  DIF_Raju$model <- input$DIF_Raju_summary_model
})
observeEvent(input$DIF_Raju_items_model, {
  DIF_Raju$model <- input$DIF_Raju_items_model
})
observeEvent(DIF_Raju$model, {
  if (DIF_Raju$model != input$DIF_Raju_summary_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_Raju_summary_model",
      selected = DIF_Raju$model
    )
  }
  if (DIF_Raju$model != input$DIF_Raju_items_model) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_Raju_items_model",
      selected = DIF_Raju$model
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_Raju_summary_correction, {
  DIF_Raju$correction <- input$DIF_Raju_summary_correction
})
observeEvent(input$DIF_Raju_items_correction, {
  DIF_Raju$correction <- input$DIF_Raju_items_correction
})
observeEvent(DIF_Raju$correction, {
  if (DIF_Raju$correction != input$DIF_Raju_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_Raju_summary_correction",
      selected = DIF_Raju$correction
    )
  }
  if (DIF_Raju$correction != input$DIF_Raju_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_Raju_items_correction",
      selected = DIF_Raju$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_Raju_summary_purification, {
  DIF_Raju$purification <- input$DIF_Raju_summary_purification
})
observeEvent(input$DIF_Raju_items_purification, {
  DIF_Raju$purification <- input$DIF_Raju_items_purification
})
observeEvent(DIF_Raju$purification, {
  if (DIF_Raju$purification != input$DIF_Raju_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_Raju_summary_purification",
      value = DIF_Raju$purification
    )
  }
  if (DIF_Raju$purification != input$DIF_Raju_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_Raju_items_purification",
      value = DIF_Raju$purification
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(binary())
  updateSliderInput(
    session = session,
    inputId = "DIF_Raju_items_item",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_Raju_method <- reactive({
  group <- unlist(group(), use.names = FALSE)
  data <- data.frame(binary())

  if (input$DIF_Raju_summary_model == "3PL") {
    guess <- itemPar3PL(data)[, 3]
  }

  fit <- tryCatch(
    switch(input$DIF_Raju_summary_model,
      "1PL" = ShinyItemAnalysis:::.difRaju_edited(
        Data = data, group = group, focal.name = 1,
        model = "1PL",
        p.adjust.method = input$DIF_Raju_summary_correction,
        purify = input$DIF_Raju_summary_purification
      ),
      "2PL" = ShinyItemAnalysis:::.difRaju_edited(
        Data = data, group = group, focal.name = 1,
        model = "2PL",
        p.adjust.method = input$DIF_Raju_summary_correction,
        purify = input$DIF_Raju_summary_purification
      ),
      "3PL" = ShinyItemAnalysis:::.difRaju_edited(
        Data = data, group = group, focal.name = 1,
        model = "3PL", c = guess,
        p.adjust.method = input$DIF_Raju_summary_correction,
        purify = input$DIF_Raju_summary_purification
      )
    ),
    error = function(e) e
  )

  validate(
    need(
      class(fit) == "Raj",
      paste0("This method cannot be used on this data. Error returned: ", fit$message)
    ),
    errorClass = "validation-error"
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY TAB ####

# ** Interpretation ####
output$DIF_Raju_summary_interpretation <- renderUI({
  type <- input$DIF_Raju_items_model

  txt <- switch(type,
    "1PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(b_{iR}\\) and \\(b_{iF}\\)
                             are difficulties for the reference and the focal group for item \\(i\\). "),
    "2PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\). "),
    "3PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters  \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\).
                             Parameter \\(c_i\\) is a common guessing parameter for item \\(i\\). ")
  )
  HTML(txt)
})

# ** Equation ####
output$DIF_Raju_summary_equation <- renderUI({
  type <- input$DIF_Raju_items_model
  eqR <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{\\theta_p - b_{iR}}}
                              {1 + e^{\\theta_p - b_{iR} }}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iR}
                              \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$")
  )
  eqF <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{\\theta_p - b_{iF}}}
                              {1 + e^{\\theta_p - b_{iF}}}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iF}
                              \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$")
  )
  paste(eqR, eqF)
})

# ** Warning for missing values ####
output$DIF_Raju_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_Raju_summary_dif_items <- renderPrint({
  DIFitems <- DIF_Raju_method()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_Raju_summary_table <- reactive({
  model <- DIF_Raju_method()

  # deal with only one pval base od model specs
  pval <- if (model$p.adjust.method == "none") {
    model$p.value
  } else {
    model$adjusted.p
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1, Inf),
    symbols = c("***", "**", "*", ".", "", "?!")
  )

  blank <- character(length(model$names))

  # IRT pars
  if (model$purification) {
    m <- nrow(model$itemParFinal) / 2

    mR <- model$itemParFinal[1:m, ]
    mF <- model$itemParFinal[(m + 1):(2 * m), ]
  } else {
    m <- nrow(model$itemParInit) / 2

    mR <- model$itemParInit[1:m, ]
    mF <- model$itemParInit[(m + 1):(2 * m), ]

    mF <- itemRescale(mR, mF) # rescaling!
  }
  par <- rbind(mR, mF)

  wh_coef <- switch(model$model,
    "1PL" = 1,
    "2PL" = 1:2,
    "3PL" = c(1, 2)
  )
  wh_sd <- switch(model$model,
    "1PL" = 2,
    "2PL" = 3:4,
    "3PL" = 3:4
  )

  R_indices <- seq(1, m)
  F_indices <- seq(m + 1, 2 * m)

  if (model$model == "3PL") {
    guess_both <- par[R_indices, 6]
  }

  coefs_R <- par[R_indices, wh_coef]
  coefs_F <- par[F_indices, wh_coef]

  if (model$model != "1PL") {
    colnames(coefs_R) <- paste0(colnames(coefs_R), "_R")
    colnames(coefs_F) <- paste0(colnames(coefs_F), "_F")
  }

  coefs <- cbind(coefs_R, coefs_F)

  if (model$model == "1PL") {
    colnames(coefs) <- c("b_R", "b_F")
  }

  se_R <- par[R_indices, wh_sd]
  se_F <- par[F_indices, wh_sd]

  if (model$model != "1PL") {
    colnames(se_R) <- paste0(colnames(se_R), "_R")
    colnames(se_F) <- paste0(colnames(se_F), "_F")
  }

  ses <- cbind(se_R, se_F)

  if (model$model == "1PL") {
    colnames(ses) <- c("se(b)_R", "se(b)_F")
  }

  pars_zigzag <-
    cbind(coefs, ses)[, order(c(seq(ncol(coefs)), seq(ncol(ses))))]
  if (model$model == "3PL") {
    pars_zigzag <- cbind(pars_zigzag, guess_both)
  }

  colnames(pars_zigzag) <- switch(model$model,
    "1PL" = c(
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))"
    ),
    "2PL" = c(
      "\\(\\mathit{a}_{R}\\)",
      "SE(\\(\\mathit{a}_{R}\\))",
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{a}_{F}\\)",
      "SE(\\(\\mathit{a}_{F}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))"
    ),
    "3PL" = c(
      "\\(\\mathit{a}_{R}\\)",
      "SE(\\(\\mathit{a}_{R}\\))",
      "\\(\\mathit{b}_{R}\\)",
      "SE(\\(\\mathit{b}_{R}\\))",
      "\\(\\mathit{a}_{F}\\)",
      "SE(\\(\\mathit{a}_{F}\\))",
      "\\(\\mathit{b}_{F}\\)",
      "SE(\\(\\mathit{b}_{F}\\))",
      "\\(\\mathit{c}\\)"
    )
  )

  # from print.raju:
  if (model$model == "1PL") {
    if (model$purification & is.null(model$anchor.names)) {
      pars <- model$itemParFinal
    } else {
      pars <- model$itemParInit
      J <- nrow(pars) / 2
      mR <- pars[1:J, 1]
      mF <- itemRescale(pars[1:J, ], pars[(J + 1):(2 * J), ])[, 1]
      rr1 <- round(mF - mR, 4)
      rr2 <- round(-2.35 * rr1, 4)
      symb1 <- symnum(abs(rr2), c(0, 1, 1.5, Inf), symbols = c("A", "B", "C"))
      matR2 <- cbind(rr1, rr2)
      # matR2 <- noquote(cbind(format(matR2, justify = "right"), symb1))
    }

    colnames(matR2) <- c("mF-mR", "deltaRaju")
  }

  tab <- if (model$model == "1PL") {
    data.frame(
      model$RajuZ,
      pval,
      pval_symb,
      blank,
      matR2[, 1],
      matR2[, 2],
      symb1, # eff size based on delta
      blank,
      pars_zigzag
    )
  } else {
    data.frame(
      model$RajuZ,
      pval,
      pval_symb,
      blank,
      pars_zigzag
    )
  }

  colnames(tab) <- if (model$model != "1PL") {
    c(
      "Raju's Z",
      ifelse(
        model$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      "",
      "",
      colnames(pars_zigzag)
    )
  } else {
    c(
      "Raju's Z",
      ifelse(
        model$p.adjust.method == "none",
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ),
      "",
      "",
      "mF-mR",
      "deltaRaju",
      "",
      "",
      colnames(pars_zigzag)
    )
  }

  rownames(tab) <- item_names()

  tab
})

output$DIF_Raju_summary_coef <- renderTable(
  {
    DIF_Raju_summary_table()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_Raju_summary_table_note <- reactive({
  model <- DIF_Raju_method()
  res <- NULL

  res$mod <- paste0("Model: ", model$model)
  res$p_adj <- paste("P-value correction method:", switch(model$p.adjust.method,
    holm = "Holm",
    hochberg = "Hochberg",
    hommel = "Hommel",
    bonferroni = "Bonferroni",
    BH = "Benjamini-Hochberg",
    BY = "Benjamini-Yekutieli",
    none = "none"
  ))
  res$puri <- paste("Item purification:", ifelse(model$purification, "used", "unutilized"))
  res$thr <- paste0(
    "Detection threshold: ", -round(model$thr, 4),
    " and ", round(model$thr, 4), " (significance level: ",
    model$alpha, ")"
  )
  res$signed <- paste0(
    "Type of Raju's Z statistic: ",
    ifelse(model$signed, "based on signed area", "based on unsigned area")
  )
  res
})

output$DIF_Raju_summary_table_note <- renderUI({
  txt_effect <- ifelse(input$DIF_Lord_summary_model == "1PL",
    "Effect size codes: 'A': negligible effect; 'B': moderate effect; 'C': large effect </br>
                       Effect size codes: 0 'A' 1.0 'B' 1.5 'C' (for absolute values of 'deltaRaju')",
    ""
  )
  HTML(
    paste(
      "Notes:",
      DIF_Raju_summary_table_note()$mod,
      DIF_Raju_summary_table_note()$p_adj,
      DIF_Raju_summary_table_note()$puri,
      DIF_Raju_summary_table_note()$thr,
      DIF_Raju_summary_table_note()$signed,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      txt_effect,
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_Raju_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_Raju_statistics.csv"
  },
  content = function(file) {
    data <- DIF_Raju_summary_table()
    data <- if (input$DIF_Raju_summary_model == "1PL") {
      data[, -c(4, 8)]
    } else {
      data[, -4]
    }

    # some dirty patch
    colnames(data) <- colnames(data) |>
      str_remove_all("\\\\\\(\\\\mathit\\{") |>
      str_remove_all("\\}\\\\\\)") |>
      str_replace("\\}_\\{", "_") |>
      str_replace("\\\\chi", "X") |>
      str_remove_all("[\\\\]")

    names(data)[3] <- "sig. symb."

    if (input$DIF_Raju_summary_model == "1PL") {
      names(data)[6] <- "eff. size"
    }

    rownames(data) <- item_names()

    write.csv(data, file)

    txt_effect <- ifelse(input$DIF_Raju_summary_model == "1PL",
      c(
        "Effect size codes: 'A': negligible effect; 'B': moderate effect; 'C': large effect",
        "Effect size codes: 0 'A' 1.0 'B' 1.5 'C' (for absolute values of 'deltaRaju')"
      ),
      ""
    )

    write(paste(
      "Notes:",
      DIF_Raju_summary_table_note()$mod,
      DIF_Raju_summary_table_note()$p_adj,
      DIF_Raju_summary_table_note()$puri,
      DIF_Raju_summary_table_note()$thr,
      DIF_Raju_summary_table_note()$signed,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      txt_effect,
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_Raju_summary_purification_info <- renderPrint({
  model <- DIF_Raju_method()
  if (input$DIF_Raju_summary_purification & !is.null(DIF_Raju_method()$difPur)) {
    cat("The table below describes purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$convergence, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_Raju_summary_purification & is.null(DIF_Raju_method()$difPur)) {
    cat("No DIF item was detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_Raju_summary_purification_table <- reactive({
  tab <- DIF_Raju_method()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_Raju_summary_purification_table <- renderTable(
  {
    DIF_Raju_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_Raju_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_Raju_purification.csv"
  },
  content = function(file) {
    data <- DIF_Raju_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_Raju_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot ####
DIF_Raju_items_plot <- reactive({
  fitRaju <- DIF_Raju_method()
  item <- input$DIF_Raju_items_item

  g <- plotDIFirt(
    parameters = fitRaju$itemParInit, test = "Raju",
    item = item, item.name = item_names()[item]
  )[[item]]
  g
})

output$DIF_Raju_items_plot <- renderPlotly({
  g <- DIF_Raju_items_plot()
  tmp1 <- g$layers[[1]]
  tmp2 <- g$layers[[2]]
  tmp3 <- g$layers[[3]]

  g$layers[[1]] <- tmp3
  g$layers[[2]] <- tmp1
  g$layers[[3]] <- tmp2

  p <- ggplotly(g)

  p$x$data[[2]]$text <- paste0(
    "Group: Reference", "<br />",
    "Ability: ", p$x$data[[2]]$x, "<br />",
    "Probability: ", p$x$data[[2]]$y
  )
  p$x$data[[3]]$text <- paste0(
    "Group: Focal", "<br />",
    "Ability: ", p$x$data[[3]]$x, "<br />",
    "Probability: ", p$x$data[[3]]$y
  )
  p$x$data[[1]]$text <- NULL

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

output$DIF_Raju_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIFIRTRaju_", item_names()[input$DIF_Raju_items_item], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_Raju_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height, width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Interpretation ####
output$DIF_Raju_items_interpretation <- renderUI({
  type <- input$DIF_Raju_items_model

  txt <- switch(type,
    "1PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(b_{iR}\\) and \\(b_{iF}\\)
                             are difficulties for the reference and the focal group for item \\(i\\). "),
    "2PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\). "),
    "3PL" = paste("As the parameters are estimated separately for the two groups, there is one
                             equation for each group. Parameters \\(a_{iR}\\) and \\(b_{iR}\\) are discrimination
                             and difficulty for the reference group for item \\(i\\). Parameters  \\(a_{iF}\\) and \\(b_{iF}\\)
                             are discrimination and difficulty for the focal group for item \\(i\\).
                             Parameter \\(c_i\\) is a common guessing parameter for item \\(i\\). ")
  )
  HTML(txt)
})

# ** Equation ####
output$DIF_Raju_items_equation <- renderUI({
  type <- input$DIF_Raju_items_model
  eqR <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{\\theta_p - b_{iR}}}
                              {1 + e^{\\theta_p - b_{iR} }}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              \\frac{e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 0\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iR}
                              \\left(\\theta_p - b_{iR} \\right)}}
                              {1 + e^{a_{iR} \\left(\\theta_p - b_{iR} \\right)}}$$")
  )
  eqF <- switch(type,
    "1PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{\\theta_p - b_{iF}}}
                              {1 + e^{\\theta_p - b_{iF}}}$$"),
    "2PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              \\frac{e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$"),
    "3PL" = paste("$$\\mathrm{P}\\left(Y_{pi} = 1 | \\theta_p, G_p = 1\\right) =
                              c_i + \\left(1 - c_i\\right) \\cdot \\frac{e^{a_{iF}
                              \\left(\\theta_p - b_{iF} \\right)}}
                              {1 + e^{a_{iF} \\left(\\theta_p - b_{iF} \\right)}}$$")
  )
  paste(eqR, eqF)
})

# ** Table of item parameters ####
DIF_Raju_items_coef <- reactive({
  fitRaju <- DIF_Raju_method()
  m <- nrow(fitRaju$itemParInit) / 2

  if (input$DIF_Raju_items_purification) {
    mR <- fitRaju$itemParFinal[1:m, ]
    mF <- fitRaju$itemParFinal[(m + 1):(2 * m), ]
  } else {
    mR <- fitRaju$itemParInit[1:m, ]
    mF <- fitRaju$itemParInit[(m + 1):(2 * m), ]
    mF <- itemRescale(mR, mF)
  }
  par <- rbind(mR, mF)

  wh_coef <- switch(input$DIF_Raju_items_model,
    "1PL" = 1,
    "2PL" = 1:2,
    "3PL" = c(1, 2, 6)
  )
  wh_sd <- switch(input$DIF_Raju_items_model,
    "1PL" = 2,
    "2PL" = 3:4,
    "3PL" = 3:4
  )
  item <- input$DIF_Raju_items_item

  tab_coef <- c(par[c(item, m + item), wh_coef])
  tab_sd <- c(par[c(item, m + item), wh_sd])

  if (input$DIF_Raju_items_model == "3PL") {
    tab_coef <- tab_coef[-6]
  }

  if (input$DIF_Raju_items_model == "3PL") {
    tab_sd <- c(tab_sd, NA)
  }

  tab <- data.frame(tab_coef, tab_sd)
  rownames(tab) <- switch(input$DIF_Raju_items_model,
    "1PL" = c("\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)"),
    "2PL" = c("\\(\\mathit{a}_{R}\\)", "\\(\\mathit{a}_{F}\\)", "\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)"),
    "3PL" = c("\\(\\mathit{a}_{R}\\)", "\\(\\mathit{a}_{F}\\)", "\\(\\mathit{b}_{R}\\)", "\\(\\mathit{b}_{F}\\)", "\\(\\mathit{c}\\)")
  )
  colnames(tab) <- c("Estimate", "SE")
  tab
})

output$DIF_Raju_items_coef <- renderTable(
  {
    DIF_Raju_items_coef()
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * METHOD COMPARISON ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

observeEvent(input$DIF_MC_unify_setting, {
  if (input$DIF_MC_matching == "score") {
    # ** Updating matching ####
    # change DMV only if custom DMV is used
    if (input$DIF_logistic_summary_matching %in% c("uploaded", "zuploaded")) {
      updateSelectInput(session,
        inputId = "DIF_logistic_summary_matching",
        selected = "zscore"
      )
    }

    if (input$DIF_NLR_summary_matching %in% c("uploaded", "zuploaded")) {
      updateSelectInput(session,
        inputId = "DIF_NLR_summary_matching",
        selected = "zscore"
      )
    }
  }

  # ** Updating purification ####
  # no need to examine already set up values
  # (when value of reactive remains the same, no additional computation is required)
  # TODO: use mapply and don't repeat the code...
  if (input$DIF_MC_purification == "purify") {
    lapply(c(
      "DIF_DP_purification",
      "DIF_MH_summary_purification",
      "DIF_logistic_summary_purification",
      "DIF_NLR_summary_purification",
      "DIF_Lord_items_purification",
      "DIF_Raju_summary_purification",
      "DIF_SIBTEST_purification"
    ), function(x) {
      updateCheckboxInput(session,
        inputId = x,
        value = TRUE
      )
    })
  } else if (input$DIF_MC_purification == "dontpurify") {
    lapply(c(
      "DIF_DP_purification",
      "DIF_MH_summary_purification",
      "DIF_logistic_summary_purification",
      "DIF_NLR_summary_purification",
      "DIF_Lord_items_purification",
      "DIF_Raju_summary_purification",
      "DIF_SIBTEST_purification"
    ), function(x) {
      updateCheckboxInput(session,
        inputId = x,
        value = FALSE
      )
    })
  }

  # ** Updating correction ####
  # note that there is no correction method in deltaPlot
  if (input$DIF_MC_correction != "asis") {
    lapply(c(
      "DIF_MH_summary_correction",
      "DIF_logistic_summary_correction",
      "DIF_NLR_summary_correction",
      "DIF_Lord_summary_correction",
      "DIF_Raju_summary_correction",
      "DIF_SIBTEST_correction"
    ), function(x) {
      updateSelectInput(session,
        inputId = x,
        selected = input$DIF_MC_correction
      )
    })
  }
})

# output to be detected by conditionalPanel which in hand shows the button
output$DIF_MC_unify_setting_condition <- reactive({
  length(DIF_MC_same_purification()) != 0 |
    length(DIF_MC_same_correction()) != 0 | length(DIF_MC_same_matching()) != 0
})

# in order to evaluate the condition above, it have to be rendered in UI
# we don't want to render any element, but we can pseudo-render it by setting:
outputOptions(output, "DIF_MC_unify_setting_condition", suspendWhenHidden = FALSE)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** MAIN TAB ####

# ** Checking purification across methods ####
DIF_MC_same_purification <- reactive({
  if (length(unique(
    list(
      input$DIF_DP_purification,
      input$DIF_MH_summary_purification,
      input$DIF_logistic_summary_purification,
      input$DIF_NLR_summary_purification,
      input$DIF_Lord_items_purification,
      input$DIF_Raju_summary_purification,
      input$DIF_SIBTEST_purification
    )
  )) != 1) {
    HTML("Warning: Purification differs across the methods.")
  }
})

output$DIF_MC_same_purification <- renderPrint({
  DIF_MC_same_purification()
})

# ** Checking correction across methods ####
DIF_MC_same_correction <- reactive({
  if (length(unique(
    list(
      input$DIF_MH_summary_correction,
      input$DIF_logistic_summary_correction,
      input$DIF_NLR_summary_correction,
      input$DIF_Lord_summary_correction,
      input$DIF_Raju_summary_correction,
      input$DIF_SIBTEST_correction
    )
  )) != 1) {
    HTML("Warning: Correction method varies across the methods.")
  }
})

output$DIF_MC_same_correction <- renderPrint({
  DIF_MC_same_correction()
})

# ** Checking matching across methods ####
DIF_MC_same_matching <- reactive({
  if (input$DIF_logistic_summary_matching %in% c("uploaded", "zuploaded") |
    input$DIF_NLR_summary_matching %in% c("uploaded", "zuploaded")) {
    HTML("Warning: Observed score should be unified across the methods!")
  }
})

output$DIF_MC_same_matching <- renderPrint({
  DIF_MC_same_matching()
})

# ** Summary table ####
output$DIF_MC_table <- renderTable(
  {
    group <- group()

    l_methods <- list()
    l_methods[["Delta"]] <- try(DIF_DP()$DIFitems)
    l_methods[["MH"]] <- try(DIF_MH_model()$DIFitems)
    l_methods[["LR"]] <- try(DIF_logistic_model()$DIFitems)
    l_methods[["NLR"]] <- try(DIF_NLR_method()$DIFitems)
    l_methods[["LORD"]] <- try(DIF_Lord_method()$DIFitems)
    l_methods[["RAJU"]] <- try(DIF_Raju_method()$DIFitems)
    l_methods[["SIBTEST"]] <- try(DIF_SIBTEST_model()$DIFitems)
    # l_methods[['DFF']] <- try(DIF_multinomial_method()$DDFitems)

    k <- length(item_names())
    idx <- lapply(l_methods, class)
    idx <- which(unlist(idx) != "try-error")

    v <- matrix(NA, ncol = length(l_methods), nrow = k)
    v[, idx] <- 0

    # there is need to handle Delta method and DDF differently
    for (j in idx) {
      if (names(l_methods)[j] == "Delta") {
        if (all(l_methods[[j]] != "no DIF item detected")) v[as.numeric(paste(l_methods[[j]])), j] <- 1
      } else {
        # if (names(l_methods)[j] == "DDF"){
        #   if (all(l_methods[[j]] != 'No DDF item detected')) v[as.numeric(paste(l_methods[[j]])), j] <- 1
        # } else {
        if (all(l_methods[[j]] != "No DIF item detected")) v[as.numeric(paste(l_methods[[j]])), j] <- 1
        # }
      }
    }

    tab <- as.data.frame(apply(v, c(1, 2), as.integer))
    rownames(tab) <- item_names()
    colnames(tab) <- names(l_methods)

    n <- nrow(tab)
    k <- ncol(tab)

    rDIF <- rowSums(tab, na.rm = TRUE)
    cDIF <- colSums(tab, na.rm = TRUE)
    cDIF[k + 1] <- 0

    tab <- cbind(tab, as.integer(rDIF))
    tab <- rbind(tab, as.integer(cDIF))

    rownames(tab)[n + 1] <- "Total"
    colnames(tab)[k + 1] <- "Total"

    tab
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# ** Note about setting below summary table ####
output$DIF_MC_table_note <- renderUI({
  HTML(paste(
    "Notes:",
    paste(
      "<b>Delta plot</b>",
      DIF_DP_table_note()$puri,
      DIF_DP_table_note()$axes,
      DIF_DP_table_note()$thr,
      sep = "<br>"
    ),
    paste(
      "<b>Mantel-Haenszel</b>",
      DIF_MH_summary_table_note()$p_adj,
      DIF_MH_summary_table_note()$puri,
      DIF_MH_summary_table_note()$thr_rounded,
      sep = "<br>"
    ),
    paste(
      "<b>Logistic regression</b>",
      DIF_logistic_summary_table_note()$matching,
      DIF_logistic_summary_table_note()$type,
      DIF_logistic_summary_table_note()$correction,
      DIF_logistic_summary_table_note()$purification,
      DIF_logistic_summary_table_note()$threshold,
      sep = "<br>"
    ),
    paste(
      "<b>Generalized logistic</b>",
      DIF_NLR_summary_table_note()$dmv,
      DIF_NLR_summary_table_note()$mod,
      DIF_NLR_summary_table_note()$type,
      DIF_NLR_summary_table_note()$p_adj,
      DIF_NLR_summary_table_note()$puri,
      DIF_NLR_summary_table_note()$thr_rounded,
      sep = "<br>"
    ),
    paste(
      "<b>IRT Lord</b>",
      DIF_Lord_summary_table_note()$mod,
      DIF_Lord_summary_table_note()$p_adj,
      DIF_Lord_summary_table_note()$puri,
      DIF_Lord_summary_table_note()$thr,
      sep = "<br>"
    ),
    paste(
      "<b>IRT Raju</b>",
      DIF_Raju_summary_table_note()$mod,
      DIF_Raju_summary_table_note()$p_adj,
      DIF_Raju_summary_table_note()$puri,
      DIF_Raju_summary_table_note()$thr,
      DIF_Raju_summary_table_note()$signed,
      sep = "<br>"
    ),
    paste(
      "<b>SIBTEST</b>",
      DIF_SIBTEST_table_note()$type,
      DIF_SIBTEST_table_note()$p_adj,
      DIF_SIBTEST_table_note()$puri,
      DIF_SIBTEST_table_note()$thr,
      sep = "<br>"
    ),
    sep = "<p>"
  ))
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * CUMULATIVE ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

match_cumulative <- c("DIF_cumulative_summary_matching", "DIF_cumulative_items_matching")
puri_cumulative <- c("DIF_cumulative_summary_purification", "DIF_cumulative_items_purification")

# ** Updating DIF matching & disable purification if DMV present ####
observe({
  if (DIF_matching_presence()) {
    lapply(match_cumulative, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore",
          "Uploaded" = "uploaded",
          "Standardized uploaded" = "zuploaded"
        ),
        selected = "zscore"
      )
    })
  } else {
    lapply(match_cumulative, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore"
        ),
        selected = "zscore"
      )
    })
  }
})

mapply(
  function(match, puri) {
    observeEvent(input[[paste0(match)]], {
      if (input[[paste0(match)]] %in% c("uploaded", "zuploaded")) {
        updateCheckboxInput(session, paste0(puri), value = FALSE)
        shinyjs::disable(paste0(puri))
      } else {
        shinyjs::enable(paste0(puri))
      }
    })
  },
  match = match_cumulative, puri = puri_cumulative
)

DIF_cumulative <- reactiveValues(
  type = NULL,
  matching = NULL,
  parametrization = NULL,
  correction = NULL,
  purification = NULL
)

# ** Updating type ####
observeEvent(input$DIF_cumulative_summary_type, {
  DIF_cumulative$type <- input$DIF_cumulative_summary_type
})
observeEvent(input$DIF_cumulative_items_type, {
  DIF_cumulative$type <- input$DIF_cumulative_items_type
})
observeEvent(DIF_cumulative$type, {
  if (DIF_cumulative$type != input$DIF_cumulative_summary_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_cumulative_summary_type",
      selected = DIF_cumulative$type
    )
  }
  if (DIF_cumulative$type != input$DIF_cumulative_items_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_cumulative_items_type",
      selected = DIF_cumulative$type
    )
  }
})

# ** Updating matching ######
observeEvent(input$DIF_cumulative_summary_matching, {
  DIF_cumulative$matching <- input$DIF_cumulative_summary_matching
})
observeEvent(input$DIF_cumulative_items_matching, {
  DIF_cumulative$matching <- input$DIF_cumulative_items_matching
})
observeEvent(DIF_cumulative$matching, {
  if (DIF_cumulative$matching != input$DIF_cumulative_summary_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_cumulative_summary_matching",
      value = DIF_cumulative$matching
    )
  }
  if (DIF_cumulative$matching != input$DIF_cumulative_items_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_cumulative_items_matching",
      value = DIF_cumulative$matching
    )
  }
})

# ** Updating parametrization ####
observeEvent(input$DIF_cumulative_summary_parametrization, {
  DIF_cumulative$parametrization <- input$DIF_cumulative_summary_parametrization
})
observeEvent(input$DIF_cumulative_items_parametrization, {
  DIF_cumulative$parametrization <- input$DIF_cumulative_items_parametrization
})
observeEvent(DIF_cumulative$parametrization, {
  if (DIF_cumulative$parametrization != input$DIF_cumulative_summary_parametrization) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_cumulative_summary_parametrization",
      selected = DIF_cumulative$parametrization
    )
  }
  if (DIF_cumulative$parametrization != input$DIF_cumulative_items_parametrization) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_cumulative_items_parametrization",
      selected = DIF_cumulative$parametrization
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_cumulative_summary_correction, {
  DIF_cumulative$correction <- input$DIF_cumulative_summary_correction
})
observeEvent(input$DIF_cumulative_items_correction, {
  DIF_cumulative$correction <- input$DIF_cumulative_items_correction
})
observeEvent(DIF_cumulative$correction, {
  if (DIF_cumulative$correction != input$DIF_cumulative_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_cumulative_summary_correction",
      selected = DIF_cumulative$correction
    )
  }
  if (DIF_cumulative$correction != input$DIF_cumulative_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_cumulative_items_correction",
      selected = DIF_cumulative$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_cumulative_summary_purification, {
  DIF_cumulative$purification <- input$DIF_cumulative_summary_purification
})
observeEvent(input$DIF_cumulative_items_purification, {
  DIF_cumulative$purification <- input$DIF_cumulative_items_purification
})
observeEvent(DIF_cumulative$purification, {
  if (DIF_cumulative$purification != input$DIF_cumulative_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_cumulative_summary_purification",
      value = DIF_cumulative$purification
    )
  }
  if (DIF_cumulative$purification != input$DIF_cumulative_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_cumulative_items_purification",
      value = DIF_cumulative$purification
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(ordinal())
  updateSliderInput(
    session = session,
    inputId = "DIF_cumulative_items",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_cumulative_method <- reactive({
  data <- ordinal()
  group <- unlist(group(), use.names = FALSE)

  data_group <- bind_cols(data, group = group)
  group_counts <- data_group |>
    drop_na() |>
    count(group)
  sc0 <- group_counts |>
    filter(group == 0) |>
    pull(n)
  sc1 <- group_counts |>
    filter(group == 1) |>
    pull(n)

  validate(
    need(sc0 > 1L, "Not enough complete observations in the reference group."),
    need(sc1 > 1L, "Not enough complete observations in the focal group.")
  )

  type <- input$DIF_cumulative_summary_type
  puri <- input$DIF_cumulative_summary_purification
  corr <- input$DIF_cumulative_summary_correction

  if (input$DIF_cumulative_summary_matching == "uploaded") {
    match <- unlist(DIFmatching())
  } else if (input$DIF_cumulative_summary_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else if (input$DIF_cumulative_summary_matching == "zscore") {
    match <- "zscore"
  } else if (input$DIF_cumulative_summary_matching == "score") {
    match <- "score"
  }

  fit <- difNLR::difORD(
    Data = data, group = group,
    focal.name = 1, model = "cumulative", match = match,
    type = type, purify = puri, p.adjust.method = corr
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY TAB ####

# ** Equation - cumulative probability ####
DIF_cumulative_summary_equation_cumulative <- reactive({
  txt1 <- "Z_p"
  txt2 <- if (input$DIF_cumulative_summary_parametrization == "irt") {
    paste0(
      "(a_{i} + a_{i\\text{DIF}} G_p)",
      "(", txt1, " - b_{ik} - b_{ik\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0k} + \\beta_{i1} ", txt1, " + \\beta_{i2} G_p + \\beta_{i3} ", txt1, ":G_p")
  }
  txt3 <- "Z_p, G_p"

  txt <- paste0("$$\\mathrm{P}(Y_{pi} \\geq k|", txt3, ") = \\frac{e^{", txt2, "}}{1 + e^{", txt2, "}}$$")
  txt
})

output$DIF_cumulative_summary_equation_cumulative <- renderUI({
  HTML(DIF_cumulative_summary_equation_cumulative())
})

# ** Equation - category probability ####
DIF_cumulative_summary_equation_category <- reactive({
  txt1 <- "Z_p"
  txt3 <- "Z_p, G_p"

  txt <- paste0("$$\\mathrm{P}(Y_{pi} = k|", txt3, ") = \\mathrm{P}(Y_{pi} \\geq k|", txt3, ") - \\mathrm{P}(Y_{pi} \\geq k + 1|", txt3, ")$$")
  txt
})

output$DIF_cumulative_summary_equation_category <- renderUI({
  HTML(DIF_cumulative_summary_equation_category())
})

# ** Warning for missing values ####
output$DIF_cumulative_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_cumulative_summary_dif_items <- renderPrint({
  DIFitems <- DIF_cumulative_method()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Function to get item parameters for cumulative and adjacent models ####
get_tidy_coefs <- function(fit) {
  par <- vctrs::vec_rbind(!!!fit$ordPAR)
  se <- vctrs::vec_rbind(!!!fit$ordSE) |>
    rename_with(~ paste0(.x, "_se"))
  res <- cbind(par, se)[, order(c(seq(ncol(par)), seq(ncol(se))))]

  res[is.na(res)] <- 0

  attr(res, "parametrization") <- fit$parametrization
  res
}

# ** Function to style summary table for cumulative and adjacent models ####
style_coef_names <- function(coef_tab, katex = TRUE) {
  # names ending with ".1" are duplicates, meaning they are SE
  nms <- names(coef_tab) |> str_replace("\\.1$", "_se")
  se <- str_detect(nms, "_se")
  idx <- str_extract(nms, "\\d+") # get one or more digits
  if (all(is.na(idx))) {
    idx <- c(1, 1)
  }

  if (attr(coef_tab, "parametrization", exact = TRUE) == "irt") {
    # parse names
    par <- str_extract(nms, "^.") # first characters
    dif <- str_detect(nms, "DIF")

    nms <- if (katex) {
      inner <- case_when(
        !is.na(idx) & dif ~ paste0(par, "_{", idx, "_{\\mathrm{DIF}}", "}"),
        is.na(idx) & dif ~ paste0(par, "_{\\mathrm{DIF}}"),
        !is.na(idx) ~ paste0(par, "_{", idx, "}"),
        TRUE ~ par
      )
      paste0("\\(\\mathit{", inner, "}\\)")
    } else {
      paste0(par, ifelse(!is.na(idx), paste0("_", idx), ""), ifelse(dif, "_DIF", ""))
    }
  } else {
    nms <- if (katex) {
      tmp <- gsub("_se", "", nms)
      tmp <- gsub("x", "\\\\beta_1", gsub("group", "\\\\beta_2", gsub("x:group", "\\\\beta_3", tmp)))
      tmp <- c(
        paste0("\\beta_{0", idx[!is.na(idx)], "}"),
        tmp[-c(1:length(idx))]
      )
      paste0("\\(\\mathit{", tmp, "}\\)")
    } else {
      tmp <- gsub("_se", "", nms)
      tmp <- gsub("x", "beta_1", gsub("group", "beta_2", gsub("x:group", "beta_3", tmp)))
      tmp <- c(
        paste0("beta_{0", idx[!is.na(idx)], "}"),
        tmp[-c(1:length(idx))]
      )
      tmp
    }
  }
  nms[se] <- paste0("SE(", nms[se], ")")
  names(coef_tab) <- nms

  coef_tab
}

# ** Summary table - statistics ####
DIF_cumulative_summary_stats <- reactive({
  fit <- DIF_cumulative_method()

  pval <- if (fit$p.adjust.method != "none") {
    fit$pval
  } else {
    fit$adj.pval
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )

  out <- data.frame(
    `LR(X^2)` = fit$Sval,
    pval,
    `sig. symb.` = pval_symb,
    check.names = FALSE
  )
  names(out)[names(out) == "pval"] <- if (fit$p.adjust.method != "none") "adj. p-value" else "p-value"

  out
})

# ** Summary table - item parameters ####
DIF_cumulative_summary_params <- reactive({
  fit <- DIF_cumulative_method()
  coefs <- coef(fit,
    SE = TRUE, IRTpars = (input$DIF_cumulative_summary_parametrization == "irt"),
    simplify = TRUE, CI = 0
  )

  estims <- coefs[c(TRUE, FALSE), ]
  ses <- coefs[c(FALSE, TRUE), ]

  out <- cbind(estims, ses)[, order(c(seq(ncol(estims)), seq(ncol(ses))))]
  attr(out, "parametrization") <- input$DIF_cumulative_summary_parametrization
  out
})

output$DIF_cumulative_summary_coef <- renderTable(
  {
    stats <- DIF_cumulative_summary_stats()
    colnames(stats) <- c(
      "LR (\\(\\mathit{\\chi^2}\\))",
      ifelse(
        "p-value" %in% names(stats),
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ), ""
    )

    params <- DIF_cumulative_summary_params()
    params <- style_coef_names(params)

    out <- data.frame(stats, "", params, check.names = FALSE, fix.empty.names = FALSE)
    rownames(out) <- item_names()
    out
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_cumulative_summary_table_note <- reactive({
  res <- NULL
  fit <- DIF_cumulative_method()

  res$matching <- paste(
    "Observed score:",
    if (fit$match[1] == "score") {
      "total score"
    } else if (fit$match[1] == "zscore") {
      "standardized total score"
    } else {
      "standardized uploaded"
    }
  )
  res$type <- paste("DIF type tested:", switch(fit$type,
    "both" = "any DIF ",
    "udif" = "uniform DIF ",
    "nudif" = "non-uniform DIF "
  ))
  res$correction <- paste("P-value correction method:", switch(fit$p.adjust.method,
    "BH" = "Benjamini-Hochberg",
    "BY" = "Benjamini-Yekutieli",
    "bonferroni" = "Bonferroni",
    "holm" = "Holm",
    "hochberg" = "Hochberg",
    "hommel" = "Hommel",
    "none" = "none"
  ))

  res
})

output$DIF_cumulative_summary_table_note <- renderUI({
  note <- DIF_cumulative_summary_table_note()

  HTML(
    paste(
      "Notes:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_cumulative_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_cumulative_table.csv"
  },
  content = function(file) {
    stats <- DIF_cumulative_summary_stats()
    params <- DIF_cumulative_summary_params()
    params <- style_coef_names(params, katex = FALSE)

    tab <- data.frame(
      item = item_names(), stats, params,
      check.names = FALSE, fix.empty.names = FALSE
    )
    note <- DIF_cumulative_summary_table_note()

    write.csv(tab, file)
    write(paste(
      "Note:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_cumulative_summary_purification_info <- renderPrint({
  fit <- DIF_cumulative_method()
  if (input$DIF_cumulative_summary_purification & !is.null(DIF_cumulative_method()$difPur)) {
    cat("The table below describes purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step,
        while the value of '0' means that the item was not detected as DIF. The first row corresponds to the initial
        classification of the items when all items were used for the calculation of the DIF matching criterion. ")
    nrIter <- fit$nrPur
    cat(
      "In this case, the convergence was", ifelse(fit$conv.puri, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_cumulative_summary_purification & is.null(DIF_cumulative_method()$difPur)) {
    cat("No DIF item was detected whatsoever, nothing to show.")
  } else {
    cat("Item purification was not requested, nothing to show.")
  }
})

# ** Purification table ####
DIF_cumulative_summary_purification_table <- reactive({
  tab <- DIF_cumulative_method()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_cumulative_summary_purification_table <- renderTable(
  {
    DIF_cumulative_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_cumulative_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_cumulative_purification.csv"
  },
  content = function(file) {
    data <- DIF_cumulative_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_cumulative_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot - cumulative ####
DIF_cumulative_items_plot_cumulative <- reactive({
  fit <- DIF_cumulative_method()
  item <- input$DIF_cumulative_items

  g <- plot(fit, item = item, plot.type = "cumulative")[[1]] +
    theme_app() +
    ggtitle(item_names()[item]) +
    theme(
      legend.box.just = "top",
      legend.justification = c("right", "bottom"),
      legend.position = "inside",
      legend.position.inside = c(0.98, 0.02),
      legend.box = "horizontal",
      legend.margin = margin(0, 0, 0, 0, unit = "cm")
    )
  g
})

output$DIF_cumulative_items_plot_cumulative <- renderPlotly({
  g <- DIF_cumulative_items_plot_cumulative()
  p <- ggplotly(g)

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- gsub("size", "Count", text)
    text <- gsub("category", "Category", text)
    text <- gsub("probability", "Probability", text)
    text <- gsub("matching", "Z-score", text)
    p$x$data[[i]]$text <- text
  }

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- lapply(strsplit(text, split = "<br />"), unique)
    text <- unlist(lapply(text, paste, collapse = "<br />"))
    p$x$data[[i]]$text <- text
  }

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot - cumulative ####
output$DIF_cumulative_items_plot_cumulative_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_cumulative_cumulative_", item_names()[input$DIF_cumulative_items], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_cumulative_items_plot_cumulative() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height,
      width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Plot - category ####
DIF_cumulative_items_plot_category <- reactive({
  fit <- DIF_cumulative_method()
  item <- input$DIF_cumulative_items

  g <- plot(fit, item = item, plot.type = "category")[[1]] +
    theme_app() +
    ggtitle(item_names()[item]) +
    theme(
      legend.box.just = "top",
      legend.justification = c("left", "top"),
      legend.position = "inside",
      legend.position.inside = c(0.02, 0.98),
      legend.box = "horizontal",
      legend.margin = margin(0, 0, 0, 0, unit = "cm")
    )
  g
})

output$DIF_cumulative_items_plot_category <- renderPlotly({
  g <- DIF_cumulative_items_plot_category()
  p <- ggplotly(g)

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- gsub("size", "Count", text)
    text <- gsub("category", "Category", text)
    text <- gsub("probability", "Probability", text)
    text <- gsub("matching", "Z-score", text)
    p$x$data[[i]]$text <- text
  }

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- lapply(strsplit(text, split = "<br />"), unique)
    text <- unlist(lapply(text, paste, collapse = "<br />"))
    p$x$data[[i]]$text <- text
  }

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot - category ####
output$DIF_cumulative_items_plot_category_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_cumulative_category_", item_names()[input$DIF_cumulative_items], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_cumulative_items_plot_category() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height,
      width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Equation - cumulative probability ####
output$DIF_cumulative_items_equation_cumulative <- renderUI({
  HTML(DIF_cumulative_summary_equation_cumulative())
})

# ** Equation - category probability ####
output$DIF_cumulative_items_equation_category <- renderUI({
  HTML(DIF_cumulative_summary_equation_category())
})

# ** Table of coefficients ####
DIF_cumulative_items_coef <- reactive({
  tab <- DIF_cumulative_summary_params()
  tab <- style_coef_names(tab)

  item <- input$DIF_cumulative_items

  tab <- t(tab[item, ])
  tab <- data.frame(tab[!grepl("SE", rownames(tab)), ], tab[grepl("SE", rownames(tab)), ])
  colnames(tab) <- c("Estimate", "SE")

  tab
})

output$DIF_cumulative_items_coef <- renderTable(
  {
    DIF_cumulative_items_coef()
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * ADJACENT ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

match_adjacent <- c("DIF_adjacent_summary_matching", "DIF_adjacent_items_matching")
puri_adjacent <- c("DIF_adjacent_summary_purification", "DIF_adjacent_items_purification")

# ** Updating DIF matching & disable purification if DMV present ####
observe({
  if (DIF_matching_presence()) {
    lapply(match_adjacent, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore",
          "Uploaded" = "uploaded",
          "Standardized uploaded" = "zuploaded"
        ),
        selected = "zscore"
      )
    })
  } else {
    lapply(match_adjacent, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore"
        ),
        selected = "zscore"
      )
    })
  }
})

mapply(
  function(match, puri) {
    observeEvent(input[[paste0(match)]], {
      if (input[[paste0(match)]] %in% c("uploaded", "zuploaded")) {
        updateCheckboxInput(session, paste0(puri), value = FALSE)
        shinyjs::disable(paste0(puri))
      } else {
        shinyjs::enable(paste0(puri))
      }
    })
  },
  match = match_adjacent, puri = puri_adjacent
)

DIF_adjacent <- reactiveValues(
  type = NULL,
  matching = NULL,
  parametrization = NULL,
  correction = NULL,
  purification = NULL
)

# ** Updating type ####
observeEvent(input$DIF_adjacent_summary_type, {
  DIF_adjacent$type <- input$DIF_adjacent_summary_type
})
observeEvent(input$DIF_adjacent_items_type, {
  DIF_adjacent$type <- input$DIF_adjacent_items_type
})
observeEvent(DIF_adjacent$type, {
  if (DIF_adjacent$type != input$DIF_adjacent_summary_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_adjacent_summary_type",
      selected = DIF_adjacent$type
    )
  }
  if (DIF_adjacent$type != input$DIF_adjacent_items_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_adjacent_items_type",
      selected = DIF_adjacent$type
    )
  }
})

# ** Updating matching ######
observeEvent(input$DIF_adjacent_summary_matching, {
  DIF_adjacent$matching <- input$DIF_adjacent_summary_matching
})
observeEvent(input$DIF_adjacent_items_matching, {
  DIF_adjacent$matching <- input$DIF_adjacent_items_matching
})
observeEvent(DIF_adjacent$matching, {
  if (DIF_adjacent$matching != input$DIF_adjacent_summary_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_summary_matching",
      value = DIF_adjacent$matching
    )
  }
  if (DIF_adjacent$matching != input$DIF_adjacent_items_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_items_matching",
      value = DIF_adjacent$matching
    )
  }
})

# ** Updating parametrization ######
observeEvent(input$DIF_adjacent_summary_parametrization, {
  DIF_adjacent$parametrization <- input$DIF_adjacent_summary_parametrization
})
observeEvent(input$DIF_adjacent_items_parametrization, {
  DIF_adjacent$parametrization <- input$DIF_adjacent_items_parametrization
})
observeEvent(DIF_adjacent$parametrization, {
  if (DIF_adjacent$parametrization != input$DIF_adjacent_summary_parametrization) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_summary_parametrization",
      value = DIF_adjacent$parametrization
    )
  }
  if (DIF_adjacent$parametrization != input$DIF_adjacent_items_parametrization) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_items_parametrization",
      value = DIF_adjacent$parametrization
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_adjacent_summary_correction, {
  DIF_adjacent$correction <- input$DIF_adjacent_summary_correction
})
observeEvent(input$DIF_adjacent_items_correction, {
  DIF_adjacent$correction <- input$DIF_adjacent_items_correction
})
observeEvent(DIF_adjacent$correction, {
  if (DIF_adjacent$correction != input$DIF_adjacent_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_adjacent_summary_correction",
      selected = DIF_adjacent$correction
    )
  }
  if (DIF_adjacent$correction != input$DIF_adjacent_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_adjacent_items_correction",
      selected = DIF_adjacent$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_adjacent_summary_purification, {
  DIF_adjacent$purification <- input$DIF_adjacent_summary_purification
})
observeEvent(input$DIF_adjacent_items_purification, {
  DIF_adjacent$purification <- input$DIF_adjacent_items_purification
})
observeEvent(DIF_adjacent$purification, {
  if (DIF_adjacent$purification != input$DIF_adjacent_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_summary_purification",
      value = DIF_adjacent$purification
    )
  }
  if (DIF_adjacent$purification != input$DIF_adjacent_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_adjacent_items_purification",
      value = DIF_adjacent$purification
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(ordinal())
  updateSliderInput(
    session = session,
    inputId = "DIF_adjacent_items",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_adjacent_model <- reactive({
  data <- ordinal()
  group <- unlist(group(), use.names = FALSE)

  data_group <- bind_cols(data, group = group)
  group_counts <- data_group |>
    drop_na() |>
    count(group)
  sc0 <- group_counts |>
    filter(group == 0) |>
    pull(n)
  sc1 <- group_counts |>
    filter(group == 1) |>
    pull(n)

  validate(
    need(sc0 > 1L, "Not enough complete observations in the reference group."),
    need(sc1 > 1L, "Not enough complete observations in the focal group.")
  )

  type <- input$DIF_adjacent_summary_type
  puri <- input$DIF_adjacent_summary_purification
  corr <- input$DIF_adjacent_summary_correction

  if (input$DIF_adjacent_summary_matching == "uploaded") {
    match <- unlist(DIFmatching())
  } else if (input$DIF_adjacent_summary_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else if (input$DIF_adjacent_summary_matching == "zscore") {
    match <- "zscore"
  } else if (input$DIF_adjacent_summary_matching == "score") {
    match <- "score"
  }

  fit <- difNLR::difORD(
    Data = data, group = group,
    focal.name = 1, model = "adjacent", match = match,
    type = type, purify = puri, p.adjust.method = corr
  )
  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY ####

# ** Equation ####
DIF_adjacent_summary_equation <- reactive({
  txt1 <- "Z_p"
  txt2 <- if (input$DIF_adjacent_summary_parametrization == "irt") {
    paste0(
      "(a_{i} + a_{i\\text{DIF}} G_p)",
      "(", txt1, " - b_{il} - b_{il\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0l} + \\beta_{i1} ", txt1, " + \\beta_{i2} G_p + \\beta_{i3} ", txt1, ":G_p")
  }
  txt3 <- "Z_p, G_p"

  txt <- paste0("$$\\mathrm{P}(Y_{pi} = k|", txt3, ") = \\frac{e^{\\sum_{l = 0}^{k}", txt2, "}}{\\sum_{j = 0}^{K_i}e^{\\sum_{l = 0}^{j}", txt2, "}}$$")
  txt
})

output$DIF_adjacent_summary_equation <- renderUI({
  HTML(DIF_adjacent_summary_equation())
})

# ** Warning for missing values ####
output$DIF_adjacent_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DIF items detected text ####
output$DIF_adjacent_summary_dif_items <- renderPrint({
  DIFitems <- DIF_adjacent_model()$DIFitems
  if (DIFitems[1] == "No DIF item detected") {
    txt <- "No item was detected as DIF."
  } else {
    txt <- paste0("Items detected as DIF items: ", paste(item_names()[DIFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table - statistics ####
DIF_adjacent_summary_stats <- reactive({
  fit <- DIF_adjacent_model()

  pval <- if (fit$p.adjust.method != "none") {
    fit$pval
  } else {
    fit$adj.pval
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )

  out <- data.frame(
    `LR(X^2)` = fit$Sval,
    pval,
    `sig. symb.` = pval_symb,
    check.names = FALSE
  )
  names(out)[names(out) == "pval"] <- if (fit$p.adjust.method != "none") "adj. p-value" else "p-value"

  out
})

# ** Summary table - item parameters ####
DIF_adjacent_summary_params <- reactive({
  fit <- DIF_adjacent_model()
  coefs <- coef(fit,
    SE = TRUE, IRTpars = (input$DIF_adjacent_summary_parametrization == "irt"),
    simplify = TRUE, CI = 0
  )

  estims <- coefs[c(TRUE, FALSE), ]
  ses <- coefs[c(FALSE, TRUE), ]

  out <- cbind(estims, ses)[, order(c(seq(ncol(estims)), seq(ncol(ses))))]
  attr(out, "parametrization") <- input$DIF_adjacent_summary_parametrization
  out
})

output$DIF_adjacent_summary_coef <- renderTable(
  {
    stats <- DIF_adjacent_summary_stats()
    colnames(stats) <- c(
      "LR (\\(\\mathit{\\chi^2}\\))",
      ifelse(
        "p-value" %in% names(stats),
        "\\(\\mathit{p}\\)-value",
        "adj. \\(\\mathit{p}\\)-value"
      ), ""
    )

    params <- DIF_adjacent_summary_params()
    params <- style_coef_names(params)

    out <- data.frame(stats, "", params, check.names = FALSE, fix.empty.names = FALSE)
    rownames(out) <- item_names()
    out
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_adjacent_summary_table_note <- reactive({
  res <- NULL
  fit <- DIF_adjacent_model()

  res$matching <- paste(
    "Observed score:",
    if (fit$match[1] == "score") {
      "total score"
    } else if (fit$match[1] == "zscore") {
      "standardized total score"
    } else {
      "standardized uploaded"
    }
  )
  res$type <- paste("DIF type tested:", switch(fit$type,
    "both" = "any DIF ",
    "udif" = "uniform DIF ",
    "nudif" = "non-uniform DIF "
  ))
  res$correction <- paste("P-value correction method:", switch(fit$p.adjust.method,
    "BH" = "Benjamini-Hochberg",
    "BY" = "Benjamini-Yekutieli",
    "bonferroni" = "Bonferroni",
    "holm" = "Holm",
    "hochberg" = "Hochberg",
    "hommel" = "Hommel",
    "none" = "none"
  ))

  res$purification <- paste("Item purification:", ifelse(fit$purification, "used", "unutilized"))
  res
})

output$DIF_adjacent_summary_table_note <- renderUI({
  note <- DIF_adjacent_summary_table_note()

  HTML(
    paste(
      "Notes:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_adjacent_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_adjacent_table.csv"
  },
  content = function(file) {
    stats <- DIF_adjacent_summary_stats()
    params <- DIF_adjacent_summary_params()
    params <- style_coef_names(params, katex = FALSE)

    tab <- data.frame(item = item_names(), stats, params, check.names = FALSE, fix.empty.names = FALSE)
    note <- DIF_adjacent_summary_table_note()

    write.csv(tab, file)

    write(paste(
      "Note:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_adjacent_summary_purification_info <- renderPrint({
  model <- DIF_adjacent_model()
  if (input$DIF_adjacent_summary_purification & !is.null(DIF_adjacent_model()$difPur)) {
    cat("The table below describes the purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DIF in (i-1)-th step, while the
        value of '0' means that the item was not detected as DIF. The first row corresponds to the initial classification of the items when all items
        were used for calculation of the DIF matching criterion. ")
    nrIter <- model$nrPur
    cat(
      "In this case, the convergence was", ifelse(model$conv.puri, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_adjacent_summary_purification & is.null(DIF_adjacent_model()$difPur)) {
    cat("No DIF items detected whatsoever, nothing to show.")
  } else {
    cat("Item purification not requested! Nothing to show.")
  }
})

# ** Purification table ####
DIF_adjacent_summary_purification_table <- reactive({
  tab <- DIF_adjacent_model()$difPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_adjacent_summary_purification_table <- renderTable(
  {
    DIF_adjacent_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_adjacent_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_adjacent_purification.csv"
  },
  content = function(file) {
    data <- DIF_adjacent_summary_purification_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_adjacent_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot ####
DIF_adjacent_items_plot <- reactive({
  fit <- DIF_adjacent_model()
  item <- input$DIF_adjacent_items

  g <- plot(fit, item = item)[[1]] +
    theme_app() +
    ggtitle(item_names()[item]) +
    theme(
      legend.box.just = "top",
      legend.justification = c("left", "top"),
      legend.position = "inside",
      legend.position.inside = c(0.02, 0.98),
      legend.box = "horizontal",
      legend.margin = margin(0, 0, 0, 0, unit = "cm")
    )
  g
})

output$DIF_adjacent_items_plot <- renderPlotly({
  g <- DIF_adjacent_items_plot()
  p <- ggplotly(g)

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- gsub("size", "Count", text)
    text <- gsub("category", "Category", text)
    text <- gsub("probability", "Probability", text)
    text <- gsub("matching", "Z-score", text)
    p$x$data[[i]]$text <- text
  }

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- lapply(strsplit(text, split = "<br />"), unique)
    text <- unlist(lapply(text, paste, collapse = "<br />"))
    p$x$data[[i]]$text <- text
  }

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Plot download ####
output$DIF_adjacent_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_adjacent_", item_names()[input$DIF_adjacent_items], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_adjacent_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height,
      width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Equation ####
output$DIF_adjacent_items_equation <- renderUI({
  HTML(DIF_adjacent_summary_equation())
})

# ** Table of coefficients ####
DIF_adjacent_items_coef <- reactive({
  tab <- DIF_adjacent_summary_params()
  tab <- style_coef_names(tab)

  item <- input$DIF_adjacent_items

  tab <- t(tab[item, ])
  tab <- data.frame(tab[!grepl("SE", rownames(tab)), ], tab[grepl("SE", rownames(tab)), ])
  colnames(tab) <- c("Estimate", "SE")

  tab
})

output$DIF_adjacent_items_coef <- renderTable(
  {
    DIF_adjacent_items_coef()
  },
  include.rownames = TRUE,
  include.colnames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * MULTINOMIAL ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** UPDATING INPUTS ####

match_multinomial <- c("DIF_multinomial_summary_matching", "DIF_multinomial_items_matching")
puri_multinomial <- c("DIF_multinomial_summary_purification", "DIF_multinomial_items_purification")

# ** Updating DIF matching & disable purification if DMV present ####
observe({
  if (DIF_matching_presence()) {
    lapply(match_multinomial, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore",
          "Uploaded" = "uploaded",
          "Standardized uploaded" = "zuploaded"
        ),
        selected = "zscore"
      )
    })
  } else {
    lapply(match_multinomial, function(i) {
      updateSelectInput(
        session,
        paste0(i),
        choices = c(
          "Total score" = "score",
          "Standardized total score" = "zscore"
        ),
        selected = "zscore"
      )
    })
  }
})

mapply(
  function(match, puri) {
    observeEvent(input[[paste0(match)]], {
      if (input[[paste0(match)]] %in% c("uploaded", "zuploaded")) {
        updateCheckboxInput(session, paste0(puri), value = FALSE)
        shinyjs::disable(paste0(puri))
      } else {
        shinyjs::enable(paste0(puri))
      }
    })
  },
  match = match_multinomial, puri = puri_multinomial
)

DIF_multinomial <- reactiveValues(
  type = NULL,
  matching = NULL,
  parametrization = NULL,
  correction = NULL,
  purification = NULL
)

# ** Updating type ####
observeEvent(input$DIF_multinomial_summary_type, {
  DIF_multinomial$type <- input$DIF_multinomial_summary_type
})
observeEvent(input$DIF_multinomial_items_type, {
  DIF_multinomial$type <- input$DIF_multinomial_items_type
})
observeEvent(DIF_multinomial$type, {
  if (DIF_multinomial$type != input$DIF_multinomial_summary_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_multinomial_summary_type",
      selected = DIF_multinomial$type
    )
  }
  if (DIF_multinomial$type != input$DIF_multinomial_items_type) {
    updateCheckboxGroupInput(
      session = session,
      inputId = "DIF_multinomial_items_type",
      selected = DIF_multinomial$type
    )
  }
})

# ** Updating matching ######
observeEvent(input$DIF_multinomial_summary_matching, {
  DIF_multinomial$matching <- input$DIF_multinomial_summary_matching
})
observeEvent(input$DIF_multinomial_items_matching, {
  DIF_multinomial$matching <- input$DIF_multinomial_items_matching
})
observeEvent(DIF_multinomial$matching, {
  if (DIF_multinomial$matching != input$DIF_multinomial_summary_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_summary_matching",
      value = DIF_multinomial$matching
    )
  }
  if (DIF_multinomial$matching != input$DIF_multinomial_items_matching) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_items_matching",
      value = DIF_multinomial$matching
    )
  }
})

# ** Updating parametrization ######
observeEvent(input$DIF_multinomial_summary_parametrization, {
  DIF_multinomial$parametrization <- input$DIF_multinomial_summary_parametrization
})
observeEvent(input$DIF_multinomial_items_parametrization, {
  DIF_multinomial$parametrization <- input$DIF_multinomial_items_parametrization
})
observeEvent(DIF_multinomial$parametrization, {
  if (DIF_multinomial$parametrization != input$DIF_multinomial_summary_parametrization) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_summary_parametrization",
      value = DIF_multinomial$parametrization
    )
  }
  if (DIF_multinomial$parametrization != input$DIF_multinomial_items_parametrization) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_items_parametrization",
      value = DIF_multinomial$parametrization
    )
  }
})

# ** Updating correction ####
observeEvent(input$DIF_multinomial_summary_correction, {
  DIF_multinomial$correction <- input$DIF_multinomial_summary_correction
})
observeEvent(input$DIF_multinomial_items_correction, {
  DIF_multinomial$correction <- input$DIF_multinomial_items_correction
})
observeEvent(DIF_multinomial$correction, {
  if (DIF_multinomial$correction != input$DIF_multinomial_summary_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_multinomial_summary_correction",
      selected = DIF_multinomial$correction
    )
  }
  if (DIF_multinomial$correction != input$DIF_multinomial_items_correction) {
    updateSelectInput(
      session = session,
      inputId = "DIF_multinomial_items_correction",
      selected = DIF_multinomial$correction
    )
  }
})

# ** Updating purification ####
observeEvent(input$DIF_multinomial_summary_purification, {
  DIF_multinomial$purification <- input$DIF_multinomial_summary_purification
})
observeEvent(input$DIF_multinomial_items_purification, {
  DIF_multinomial$purification <- input$DIF_multinomial_items_purification
})
observeEvent(DIF_multinomial$purification, {
  if (DIF_multinomial$purification != input$DIF_multinomial_summary_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_summary_purification",
      value = DIF_multinomial$purification
    )
  }
  if (DIF_multinomial$purification != input$DIF_multinomial_items_purification) {
    updateCheckboxInput(
      session = session,
      inputId = "DIF_multinomial_items_purification",
      value = DIF_multinomial$purification
    )
  }
})

# ** Updating item slider ####
observe({
  item_count <- ncol(ordinal())
  updateSliderInput(
    session = session,
    inputId = "DIF_multinomial_items",
    max = item_count
  )
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** METHOD ####

# ** Method fit ####
DIF_multinomial_method <- reactive({
  data <- nominal()
  group <- unlist(group())
  key <- key()

  type <- input$DIF_multinomial_summary_type
  puri <- input$DIF_multinomial_summary_purification
  corr <- input$DIF_multinomial_summary_correction

  if (input$DIF_multinomial_summary_matching == "uploaded") {
    match <- unlist(DIFmatching())
  } else if (input$DIF_multinomial_summary_matching == "zuploaded") {
    match <- scale(apply(as.data.frame(unlist(DIFmatching())), 1, sum))
  } else if (input$DIF_multinomial_summary_matching == "zscore") {
    match <- "zscore"
  } else if (input$DIF_multinomial_summary_matching == "score") {
    match <- "score"
  }

  fit <- tryCatch(
    ddfMLR(
      Data = data, group = group, focal.name = 1, match = match,
      key = key, p.adjust.method = corr,
      type = type, purify = puri
    ),
    error = function(e) e
  )

  validate(
    need(
      class(fit) == "ddfMLR",
      paste0("This method cannot be used on this data. Error returned: ", fit$message)
    ),
    errorClass = "validation-error"
  )

  fit
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** SUMMARY ####

# ** Equation - correct answer ####
DIF_multinomial_summary_equation_correct <- reactive({
  txt1 <- "Z_p"
  txt2 <- if (input$DIF_multinomial_summary_parametrization == "irt") {
    paste0(
      "(a_{il} + a_{il\\text{DIF}} G_p)",
      "(", txt1, " - b_{il} - b_{il\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0l} + \\beta_{i1l} ", txt1, " + \\beta_{i2l} G_p + \\beta_{i3l} ", txt1, ":G_p")
  }
  txt3 <- "Z_p, G_p"

  txt <- paste0("$$\\mathrm{P}(Y_{pi} = K_i|", txt3, ") = \\frac{1}{\\sum_{l} e^{", txt2, "}}$$")
  txt
})

output$DIF_multinomial_summary_equation_correct <- renderUI({
  HTML(DIF_multinomial_summary_equation_correct())
})

# ** Equation - distractor ####
DIF_multinomial_summary_equation_distractor <- reactive({
  txt1 <- "Z_p"
  txt2 <- if (input$DIF_multinomial_summary_parametrization == "irt") {
    paste0(
      "(a_{ik} + a_{ik\\text{DIF}} G_p)",
      "(", txt1, " - b_{ik} - b_{ik\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0k} + \\beta_{i1k} ", txt1, " + \\beta_{i2k} G_p + \\beta_{i3k} ", txt1, ":G_p")
  }
  txt3 <- if (input$DIF_multinomial_summary_parametrization == "irt") {
    paste0(
      "(a_{il} + a_{il\\text{DIF}} G_p)",
      "(", txt1, " - b_{il} - b_{il\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0l} + \\beta_{i1l} ", txt1, " + \\beta_{i2l} G_p + \\beta_{i3l} ", txt1, ":G_p")
  }
  txt4 <- "Z_p, G_p"

  txt <- paste0("$$\\mathrm{P}(Y_{pi} = k|", txt4, ") = \\frac{e^{", txt2, "}}{\\sum_{l} e^{", txt3, "}}$$")
  txt
})

output$DIF_multinomial_summary_equation_distractor <- renderUI({
  HTML(DIF_multinomial_summary_equation_distractor())
})

# ** Warning for missing values ####
output$DIF_multinomial_summary_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** DDF items detected text ####
output$DIF_multinomial_summary_dif_items <- renderPrint({
  DDFitems <- DIF_multinomial_method()$DDFitems
  if (DDFitems[1] == "No DDF item detected") {
    txt <- "No item was detected as DDF."
  } else {
    txt <- paste0("Items detected as DDF items: ", paste(item_names()[DDFitems], collapse = ", "))
  }
  HTML(txt)
})

# ** Summary table ####
DIF_multinomial_summary_coef <- reactive({
  fit <- DIF_multinomial_method()

  # only one pval var, based on model specs
  pval <- if (fit$p.adjust.method == "none") {
    fit$pval
  } else {
    fit$adj.pval
  }

  pval_symb <- symnum(pval,
    c(0, 0.001, 0.01, 0.05, 0.1, 1),
    symbols = c("***", "**", "*", ".", "")
  )
  blank <- character(ncol(fit$Data))

  # table
  tab_stat <-
    data.frame(
      check.names = FALSE,
      fit$Sval,
      pval,
      pval_symb
    )

  colnames(tab_stat) <- c(
    "LR (\\(\\mathit{\\chi^2}\\))",
    ifelse(
      fit$p.adjust.method == "none",
      "\\(\\mathit{p}\\)-value",
      "adj. \\(\\mathit{p}\\)-value"
    ),
    ""
  )

  rownames(tab_stat) <- item_names()
  attr(tab_stat, "pval_adj") <- isTRUE(fit$p.adjust.method != "none")
  tab_stat
})

output$DIF_multinomial_summary_coef <- renderTable(
  {
    DIF_multinomial_summary_coef()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Note about setting below summary table ####
DIF_multinomial_summary_table_note <- reactive({
  res <- NULL
  fit <- DIF_multinomial_method()

  res$matching <- paste(
    "DDF matching variable:",
    if (fit$match[1] == "score") {
      "total score"
    } else if (fit$match[1] == "zscore") {
      "standardized total score"
    } else {
      "standardized uploaded"
    }
  )
  res$type <- paste("DDF type tested:", switch(fit$type,
    "both" = "any DDF ",
    "udif" = "uniform DDF ",
    "nudif" = "non-uniform DDF "
  ))
  res$correction <- paste("P-value correction method:", switch(fit$p.adjust.method,
    "BH" = "Benjamini-Hochberg",
    "BY" = "Benjamini-Yekutieli",
    "bonferroni" = "Bonferroni",
    "holm" = "Holm",
    "hochberg" = "Hochberg",
    "hommel" = "Hommel",
    "none" = "none"
  ))
  res$purification <- paste("Item purification:", ifelse(fit$purification, "used", "unutilized"))
  res
})

output$DIF_multinomial_summary_table_note <- renderUI({
  note <- DIF_multinomial_summary_table_note()

  HTML(
    paste(
      "Notes:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "</br>"
    )
  )
})

# ** Download summary table ####
output$DIF_multinomial_summary_table_download <- downloadHandler(
  filename = function() {
    "DIF_multinomial_table.csv"
  },
  content = function(file) {
    table <- DIF_multinomial_summary_coef()

    colnames(table) <- c(
      "LR (X^2})",
      ifelse(
        attr(table, "pval_adj"),
        "adj. p-value",
        "p-value"
      ),
      "sig. symb."
    )

    note <- DIF_multinomial_summary_table_note()

    write.csv(table, file)
    write(paste(
      "Note:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Item parameters table ####
DIF_multinomial_summary_coef_parameters <- reactive({
  fit <- DIF_multinomial_method()

  # estimated coefficients for all items with standard errors
  coefs <- coef(fit, SE = TRUE, IRTpars = (input$DIF_multinomial_summary_parametrization == "irt"), simplify = TRUE, CI = 0)

  # adding missing columns if necessary / ordering columns
  if (input$DIF_multinomial_summary_parametrization == "classic") {
    new_cols <- c(group = 0, "x:group" = 0)
    coefs <- coefs |> add_column(!!!new_cols[setdiff(names(new_cols), names(coefs))])
  } else {
    coefs <- coefs[, c(which(grepl("a", colnames(coefs))), which(!grepl("a", colnames(coefs))))]
  }
  estims <- coefs[grepl("estimate", rownames(coefs)), ]
  ses <- coefs[!grepl("estimate", rownames(coefs)), ]

  pars <- cbind(estims, ses)[, order(c(seq(ncol(estims)), seq(ncol(ses))))]

  if (input$DIF_multinomial_summary_parametrization == "irt") {
    est_pars <- c("a", "a_{\\mathrm{DIF}}", "b", "b_{\\mathrm{DIF}}")
    n_est_pars <- ncol(pars) / 2
    which_est_pars <- switch(paste(n_est_pars),
      "2" = c(1, 3),
      "3" = c(1, 3, 4),
      "4" = 1:4
    )
    est_pars <- est_pars[which_est_pars]

    colnames(pars) <- paste0(
      c("", "SE("),
      paste0("\\(\\mathit{", rep(est_pars, each = 2), "}\\)"),
      c("", ")")
    )
  } else {
    est_pars <- c("\\beta_0", "\\beta_1", "\\beta_2", "\\beta_3")
    n_est_pars <- ncol(pars) / 2
    which_est_pars <- switch(paste(n_est_pars),
      "2" = 1:2,
      "3" = 1:3,
      "4" = 1:4
    )
    est_pars <- est_pars[which_est_pars]

    colnames(pars) <- paste0(
      c("", "SE("),
      paste0("\\(\\mathit{", rep(est_pars, each = 2), "}\\)"),
      c("", ")")
    )
  }

  rownames(pars) <- gsub("estimate_", "", rownames(pars))
  pars
})

output$DIF_multinomial_summary_coef_parameters <- renderTable(
  {
    DIF_multinomial_summary_coef_parameters()
  },
  rownames = TRUE,
  colnames = TRUE
)

# ** Download item parameters table ####
output$DIF_multinomial_summary_parameters_download <- downloadHandler(
  filename = function() {
    "DIF_multinomial_parameters.csv"
  },
  content = function(file) {
    table <- DIF_multinomial_summary_coef_parameters()

    # remove math
    names(table) <- names(table) |>
      str_remove_all("\\\\\\(\\\\mathit\\{") |>
      str_remove_all("\\\\mathrm\\{") |>
      str_remove_all("\\\\\\)") |>
      str_remove_all("[\\\\{}}]")

    note <- DIF_multinomial_summary_table_note()

    write.csv(table, file)
    write(paste(
      "Note:",
      note$matching,
      note$type,
      note$correction,
      note$purification,
      "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
      sep = "\n"
    ), file, append = TRUE)
  }
)

# ** Purification info - number of iterations ####
output$DIF_multinomial_summary_purification_info <- renderPrint({
  fit <- DIF_multinomial_method()
  if (input$DIF_multinomial_summary_purification & !is.null(DIF_multinomial_method()$ddfPur)) {
    cat("The table below describes the purification process. The rows correspond to the purification iteration and the columns
        correspond to items. Value of '1' in the i-th row means that an item was detected as DDF in (i-1)-th step, while the
        value of '0' means that the item was not detected as DDF. The first row corresponds to the initial classification of the items when all items
        were used for calculation of the DIF matching criterion. ")
    nrIter <- fit$nrPur
    cat(
      "In this case, the convergence was", ifelse(fit$conv.puri, "reached", "NOT reached even"), "after", nrIter,
      ifelse(nrIter == 1, "iteration.", "iterations.")
    )
  } else if (input$DIF_multinomial_summary_purification & is.null(DIF_multinomial_method()$ddfPur)) {
    cat("No DDF items detected whatsoever, nothing to show.")
  } else {
    cat("Item purification not requested! Nothing to show.")
  }
})

# ** Purification table ####
DIF_multinomial_summary_purification_table <- reactive({
  tab <- DIF_multinomial_method()$ddfPur

  if (!is.null(tab)) {
    colnames(tab) <- item_names()
    rownames(tab) <- paste0("Step ", seq(0, nrow(tab) - 1))
    tab
  }
})

output$DIF_multinomial_summary_purification_table <- renderTable(
  {
    DIF_multinomial_summary_purification_table()
  },
  rownames = TRUE,
  colnames = TRUE,
  digits = 0
)

# ** Download purification table ####
output$DIF_multinomial_summary_purification_table_download <- downloadHandler(
  filename = function() {
    "DIF_multinomial_purification.csv"
  },
  content = function(file) {
    data <- dif_multinomial_puri_table()
    write.csv(data, file)
  }
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** ITEMS ####

# ** Warning for missing values ####
output$DIF_multinomial_items_na_alert <- renderUI({
  txt <- na_score()
  HTML(txt)
})

# ** Plot all items ####
DIF_multinomial_items_plot_all <- reactive({
  fit <- DIF_multinomial_method()

  g <- list()
  for (i in 1:ncol(nominal())) {
    g[[i]] <- plot(fit, item = i)[[1]] +
      theme_app() +
      ggtitle(item_names()[i]) +
      theme(
        legend.box.just = "top",
        legend.justification = c("left", "top"),
        legend.position = "inside",
        legend.position.inside = c(0.02, 0.98),
        legend.box = "horizontal",
        legend.margin = margin(0, 0, 0, 0, unit = "cm")
      )
  }
  g
})

# ** Plot ####
DIF_multinomial_items_plot <- reactive({
  item <- input$DIF_multinomial_items

  g <- DIF_multinomial_items_plot_all()
  g[[item]]
})

output$DIF_multinomial_items_plot <- renderPlotly({
  g <- DIF_multinomial_items_plot()
  p <- ggplotly(g)

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- gsub("answ", "Category", text)
    text <- gsub("variable", "Category", text)
    text <- gsub("Freq.1", "Count", text)
    text <- gsub("Freq", "Empirical probability", text)
    text <- gsub("value", "Probability", text)
    text <- gsub("score", "Z-score", text)
    p$x$data[[i]]$text <- text
  }

  for (i in 1:length(p$x$data)) {
    text <- p$x$data[[i]]$text
    text <- lapply(strsplit(text, split = "<br />"), unique)
    text <- unlist(lapply(text, paste, collapse = "<br />"))
    p$x$data[[i]]$text <- text
  }

  p$elementId <- NULL
  hide_legend(p) |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot ####
output$DIF_multinomial_items_plot_download <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_multinomial_", item_names()[input$DIF_multinomial_items], ".png")
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_multinomial_items_plot() +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height,
      width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Download all plots ####
output$DIF_multinomial_items_plot_download_all <- downloadHandler(
  filename = function() {
    paste0("fig_DIF_multinomial_all.zip")
  },
  content = function(file) {
    # go to a temp dir to avoid permission issues
    owd <- setwd(tempdir())
    on.exit(setwd(owd))
    files <- NULL

    # loop through the sheets
    for (i in 1:ncol(ordinal())) {
      # write each sheet to a csv file, save the name
      fileName <- paste0("fig_DIF_multinomial_", item_names()[i], ".png")
      ggsave(fileName,
        plot = DIF_multinomial_items_plot_all()[[i]] +
          theme(text = element_text(size = setting_figures$text_size)),
        device = "png",
        height = setting_figures$height,
        width = setting_figures$width,
        dpi = setting_figures$dpi
      )
      files <- c(fileName, files)
    }
    # create the zip file
    zip(file, files)
  }
)

# ** Equation ####
output$DIF_multinomial_items_equation <- renderUI({
  item <- input$DIF_multinomial_items
  key <- key()

  txt1 <- "Z_p"
  txt2 <- if (input$DIF_multinomial_items_parametrization == "irt") {
    paste0(
      "(a_{ik} + a_{ik\\text{DIF}} G_p)",
      "(", txt1, " - b_{ik} - b_{ik\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0k} + \\beta_{i1k} ", txt1, " + \\beta_{i2k} G_p + \\beta_{i3k} ", txt1, ":G_p")
  }
  txt3 <- if (input$DIF_multinomial_summary_parametrization == "irt") {
    paste0(
      "(a_{il} + a_{il\\text{DIF}} G_p)",
      "(", txt1, " - b_{il} - b_{il\\text{DIF}} G_p)"
    )
  } else {
    paste0("\\beta_{i0l} + \\beta_{i1l} ", txt1, " + \\beta_{i2l} G_p + \\beta_{i3l} ", txt1, ":G_p")
  }
  txt4 <- "Z_p, G_p"


  cor_option <- key[item]
  HTML(paste(sprintf(
    paste("For item %s, corresponding equations for the multinomial model are given by:
      $$\\mathrm{P}(Y_{pi} = %s|", txt4, ") = \\frac{1}{\\sum_{l} e^{", txt3, "}}, $$,
      $$\\mathrm{P}(Y_{pi} = k|", txt4, ") = \\frac{e^{", txt2, "}}{\\sum_{l} e^{", txt3, "}}, $$
      where \\(%s\\) is the correct answer and \\(k\\) is one of the wrong options (distractors). "),
    item, cor_option, cor_option
  )))
})

# ** Table of coefficients ####
output$DIF_multinomial_items_coef <- renderTable(
  {
    item <- input$DIF_multinomial_items
    item_names <- gsub(" ", "\\.", item_names())

    tab <- DIF_multinomial_summary_coef_parameters()
    tab <- tab[grepl(paste(item_names[item], ""), rownames(tab)), ]

    rownames(tab) <- gsub(item_names[item], "", rownames(tab))
    tab
  },
  rownames = TRUE
)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# ** REPORTS ####

# ** Method ####
report_DIF_multinomial_method <- reactive({
  if (!input$customizeCheck) {
    fit <- DIF_multinomial_method()
  } else {
    data <- nominal()
    group <- unlist(group())
    key <- key()

    type <- input$type_DIF_report
    puri <- input$puri_DIF_report
    corr <- input$correction_method_DIF_report

    fit <- tryCatch(
      ddfMLR(
        Data = data, group = group, focal.name = 1, # match = match,
        key = key, p.adjust.method = corr,
        type = type, purify = puri, parametrization = "classic"
      ),
      error = function(e) e
    )

    validate(
      need(
        class(fit) == "ddfMLR",
        paste0("This method cannot be used on this data. Error returned: ", fit$message)
      ),
      errorClass = "validation-error"
    )
  }

  fit
})

# ** Plot ####
report_DIF_multinomial_plot <- reactive({
  fit <- report_DIF_multinomial_method()

  graflist <- list()
  if (fit$DDFitems[[1]] != "No DDF item detected") {
    for (i in 1:length(fit$DDFitems)) {
      g <- plot(fit, item = fit$DDFitems[i])[[1]] +
        theme(
          text = element_text(size = 12),
          plot.title = element_text(size = 12, face = "bold", vjust = 1.5)
        ) +
        ggtitle(paste("\nDDF multinomial plot for item ", item_numbers()[fit$DDFitems[i]]))
      graflist[[i]] <- g
    }
  } else {
    graflist <- NULL
  }
  graflist
})

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * TRAINING ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ####
# * MAIN TAB ####

# *** Interpretation ####
output$DIF_training_interpretation <- renderUI({
  aR <- input$DIF_training_parameter_aR
  bR <- input$DIF_training_parameter_bR
  cR <- 0
  dR <- 1

  aF <- input$DIF_training_parameter_aF
  bF <- input$DIF_training_parameter_bF
  cF <- 0
  dF <- 1

  theta0 <- input$DIF_training_parameter_theta

  ccirt <- function(theta, a, b, c, d) {
    return(c + (d - c) / (1 + exp(-a * (theta - b))))
  }

  probR <- ccirt(theta0, a = aR, b = bR, c = cR, d = dR)
  probF <- ccirt(theta0, a = aF, b = bF, c = cF, d = dF)

  txt1 <- paste0(
    "In the <font color='blue'><b>reference</b></font> group, a respondent with
    the ability ",
    paste0("\\(\\theta= ", theta0, "\\)"),
    " has the probability of the correct answer to an item with parameters ",
    paste0("\\(a = ", aR, "\\),"), " ",
    paste0("\\(b = ", bR, "\\),"), " ",
    paste0("\\(c = ", cR, "\\),"), " and ",
    paste0("\\(d = ", dR, "\\)"),
    " equal to <b>", sprintf("%.2f", probR), "</b>. "
  )
  txt2 <- paste0(
    "In the <font color='#e6b800'><b>focal</b></font> group, a respondent with
    the ability ",
    paste0("\\(\\theta= ", theta0, "\\)"),
    "has the probability of the correct answer  to an item with parameters ",
    paste0("\\(a = ", aF, "\\),"), " ",
    paste0("\\(b = ", bF, "\\),"), " ",
    paste0("\\(c = ", cF, "\\),"), " and ",
    paste0("\\(d = ", dF, "\\)"),
    " equal to <b>", sprintf("%.2f", probF), "</b>. "
  )

  txt <- paste0("<b>Interpretation: </b>", txt1, txt2)
  HTML(txt)
})

# ** Plot for training ####
DIF_training_plot <- reactive({
  aR <- input$DIF_training_parameter_aR
  bR <- input$DIF_training_parameter_bR
  cR <- 0
  dR <- 1

  aF <- input$DIF_training_parameter_aF
  bF <- input$DIF_training_parameter_bF
  cF <- 0
  dF <- 1

  theta0 <- input$DIF_training_parameter_theta

  ccirt <- function(theta, a, b, c, d) {
    return(c + (d - c) / (1 + exp(-a * (theta - b))))
  }

  probR <- ccirt(theta0, a = aR, b = bR, c = cR, d = dR)
  probF <- ccirt(theta0, a = aF, b = bF, c = cF, d = dF)

  df <- data.frame(
    Reference = ccirt(seq(-4, 4, 0.01), aR, bR, cR, dR),
    Focal = ccirt(seq(-4, 4, 0.01), aF, bF, cF, dF),
    Ability = seq(-4, 4, 0.01)
  )

  df <- tidyr::pivot_longer(df, -Ability, names_to = "Group", values_to = "Probability")
  df$Group <- factor(df$Group, levels = c("Reference", "Focal"))

  g <- ggplot(data = df, aes(x = Ability, y = Probability, col = Group, linetype = Group)) +
    geom_line(size = 0.8) +
    geom_segment(
      aes(
        y = probR, yend = probR,
        x = -4, xend = theta0
      ),
      color = "gray", linetype = "dotdash"
    ) +
    geom_segment(
      aes(
        y = probF, yend = probF,
        x = -4, xend = theta0
      ),
      color = "gray", linetype = "dotdash"
    ) +
    geom_segment(
      aes(
        y = 0,
        yend = max(probR, probF),
        x = theta0, xend = theta0
      ),
      color = "gray", linetype = "dotdash"
    ) +
    xlim(-4, 4) +
    xlab("Ability") +
    ylab("Probability of correct answer") +
    ylim(0, 1) +
    scale_color_manual(
      name = "",
      values = c("blue", "#e6b800"),
      labels = c(
        paste(paste(letters[1:2], "=", c(aR, bR)),
          collapse = ", "
        ),
        paste(paste(paste(letters[1:2], "=", c(aF, bF))),
          collapse = ", "
        )
      )
    ) +
    scale_linetype_manual(
      name = "",
      values = c("solid", "dashed"),
      labels = c(
        paste(paste(letters[1:2], "=", c(aR, bR)),
          collapse = ", "
        ),
        paste(paste(paste(letters[1:2], "=", c(aF, bF))),
          collapse = ", "
        )
      )
    ) +
    theme_app()
  g
})

output$DIF_training_plot <- renderPlotly({
  g <- DIF_training_plot()
  p <- ggplotly(g)

  # Reference group, probabilities
  text <- gsub("Group: Reference<br />Group: Reference", "Group: Reference", p$x$data[[1]]$text)
  p$x$data[[1]]$text <- text

  # Focal group, probabilities
  text <- gsub("Group: Focal<br />Group: Focal", "Group: Focal", p$x$data[[2]]$text)
  p$x$data[[2]]$text <- text

  # Reference group, selected theta
  text <- paste(strsplit(p$x$data[[3]]$text, "<br />")[[1]][-c(1:4)], collapse = "<br />")
  text <- gsub("Group: gray<br />Group: dotdash", "Group: Reference", text)
  p$x$data[[3]]$text <- text

  # Focal group, selected theta
  text <- paste(strsplit(p$x$data[[4]]$text, "<br />")[[1]][-c(1:4)], collapse = "<br />")
  text <- gsub("Group: gray<br />Group: dotdash", "Group: Focal", text)
  p$x$data[[4]]$text <- text

  # Selected theta
  text <- paste(strsplit(p$x$data[[5]]$text, "<br />")[[1]][5], collapse = "<br />")
  text <- paste(text, "<br />Group: Reference<br />Group: Focal")
  p$x$data[[5]]$text <- text

  p$elementId <- NULL

  p |> plotly::config(displayModeBar = FALSE)
})

# ** Download plot ####
output$DIF_training_plot_download <- downloadHandler(
  filename = function() {
    "fig_DIFtraining.png"
  },
  content = function(file) {
    ggsave(file,
      plot = DIF_training_plot() +
        theme(
          legend.position = "inside",
          legend.position.inside = c(0.97, 0.03),
          legend.justification = c(0.97, 0.03)
        ) +
        theme(text = element_text(size = setting_figures$text_size)),
      device = "png",
      height = setting_figures$height,
      width = setting_figures$width,
      dpi = setting_figures$dpi
    )
  }
)

# ** Exercise 1 ####
# *** Correct answers for Exercise 1 ####
DIF_training_correct_answers_1 <- reactive({
  ccirt <- function(theta, a, b, c, d) {
    return(1 / (1 + exp(-a * (theta - b))))
  }
  # Exercise 1.1
  # correct parameters for reference and focal group
  aR <- 1
  bR <- 0
  aF <- 1
  bF <- 1
  parR <- c(aR, bR)
  parF <- c(aF, bF)

  # Exercise 1.3
  # probability calculation
  theta0 <- c(-2, 0, 2)

  probR <- ccirt(theta0, aR, bR)
  probF <- ccirt(theta0, aF, bF)

  correct_answers <- list(
    Ex1_1 = list(
      parR = parR,
      parF = parF
    ),
    Ex1_2 = "uniform",
    Ex1_3 = list(
      probR = probR,
      probF = probF
    ),
    Ex1_4 = "reference"
  )
  correct_answers
})

# ** Evaluation of answers for Exercise 1 ####
DIF_training_answers_check_1 <- eventReactive(input$DIF_training_1_submit, {
  correct_answers <- DIF_training_correct_answers_1()

  # Exercise 1.1
  aR <- input$DIF_training_parameter_aR
  bR <- input$DIF_training_parameter_bR
  aF <- input$DIF_training_parameter_aF
  bF <- input$DIF_training_parameter_bF
  parR <- c(aR, bR)
  parF <- c(aF, bF)

  check1_1 <- c(all(abs(parR - correct_answers[["Ex1_1"]]$parR) <= 0.05) &
    all(abs(parF - correct_answers[["Ex1_1"]]$parF) <= 0.05))
  # Exercise 1.2
  check1_2 <- input$DIF_training_1_2 == correct_answers[["Ex1_2"]]

  # Exercise 1.3
  probR <- c(input$DIF_training_1_3_1R, input$DIF_training_1_3_2R, input$DIF_training_1_3_3R)
  probF <- c(input$DIF_training_1_3_1F, input$DIF_training_1_3_2F, input$DIF_training_1_3_3F)
  check1_3R <- abs(probR - correct_answers[["Ex1_3"]]$probR) <= 0.05
  check1_3F <- abs(probF - correct_answers[["Ex1_3"]]$probF) <= 0.05

  # Exercise 1.4
  check1_4 <- input$DIF_training_1_4 == correct_answers[["Ex1_4"]]

  check <- list(
    check1_1 = check1_1,
    check1_2 = check1_2,
    check1_3R = check1_3R,
    check1_3F = check1_3F,
    check1_4 = check1_4
  )
  res <- sum(sapply(check, sum)) / sum(sapply(check, length))
  ans <- lapply(check, function(x) {
    ifelse(is.na(x),
      "<b><font color = 'red'>!</font></b>",
      ifelse(x,
        "<font color='green'>&#10004;</font>",
        "<font color='red'>&#10006;</font>"
      )
    )
  })
  ans[["total"]] <- res
  ans
})

# *** Checkmarks for Exercise 1 ####
output$DIF_training_1_1_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_1"]])
})
output$DIF_training_1_2_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_2"]])
})
output$DIF_training_1_3_1R_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3R"]][1])
})
output$DIF_training_1_3_2R_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3R"]][2])
})
output$DIF_training_1_3_3R_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3R"]][3])
})
output$DIF_training_1_3_1F_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3F"]][1])
})
output$DIF_training_1_3_2F_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3F"]][2])
})
output$DIF_training_1_3_3F_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_3F"]][3])
})
output$DIF_training_1_4_answer <- renderUI({
  HTML(DIF_training_answers_check_1()[["check1_4"]])
})
output$DIF_training_1_answer <- renderUI({
  res <- DIF_training_answers_check_1()[["total"]]
  HTML(ifelse(is.na(res),
    "<font color = 'red'>Check the format</font>",
    ifelse(res == 1,
      "<font color='green'>Everything correct! Well done!</font>",
      paste0("<font color='red'>", round(100 * res), "% correct. Try again.</font>")
    )
  ))
})

# ** Exercise 2 ####
# *** Correct answers for Exercise 2 ####
DIF_training_correct_answers_2 <- reactive({
  ccirt <- function(theta, a, b, c, d) {
    return(1 / (1 + exp(-a * (theta - b))))
  }
  # Exercise 2.1
  # correct parameters for reference and focal group
  aR <- 0.8
  bR <- -0.5
  aF <- 1.5
  bF <- 1
  parR <- c(aR, bR)
  parF <- c(aF, bF)

  # Exercise 2.3
  # probability calculation
  theta0 <- c(-1, 0, 1)

  probR <- ccirt(theta0, aR, bR)
  probF <- ccirt(theta0, aF, bF)

  correct_answers <- list(
    Ex2_1 = list(
      parR = parR,
      parF = parF
    ),
    Ex2_2 = "nonuniform",
    Ex2_3 = list(
      probR = probR,
      probF = probF
    ),
    Ex2_4 = "depends"
  )
  correct_answers
})

# ** Evaluation of answers for Exercise 2 ####
DIF_training_answers_check_2 <- eventReactive(input$DIF_training_2_submit, {
  correct_answers <- DIF_training_correct_answers_2()

  # Exercise 2.1
  aR <- input$DIF_training_parameter_aR
  bR <- input$DIF_training_parameter_bR
  aF <- input$DIF_training_parameter_aF
  bF <- input$DIF_training_parameter_bF
  parR <- c(aR, bR)
  parF <- c(aF, bF)

  check2_1 <- c(all(abs(parR - correct_answers[["Ex2_1"]]$parR) <= 0.05) &
    all(abs(parF - correct_answers[["Ex2_1"]]$parF) <= 0.05))
  # Exercise 2.2
  check2_2 <- input$DIF_training_2_2 == correct_answers[["Ex2_2"]]

  # Exercise 2.3
  probR <- c(input$DIF_training_2_3_1R, input$DIF_training_2_3_2R, input$DIF_training_2_3_3R)
  probF <- c(input$DIF_training_2_3_1F, input$DIF_training_2_3_2F, input$DIF_training_2_3_3F)
  check2_3R <- abs(probR - correct_answers[["Ex2_3"]]$probR) <= 0.05
  check2_3F <- abs(probF - correct_answers[["Ex2_3"]]$probF) <= 0.05

  # Exercise 2.4
  check2_4 <- input$DIF_training_2_4 == correct_answers[["Ex2_4"]]

  check <- list(
    check2_1 = check2_1,
    check2_2 = check2_2,
    check2_3R = check2_3R,
    check2_3F = check2_3F,
    check2_4 = check2_4
  )
  res <- sum(sapply(check, sum)) / sum(sapply(check, length))
  ans <- lapply(check, function(x) {
    ifelse(is.na(x),
      "<b><font color = 'red'>!</font></b>",
      ifelse(x,
        "<font color='green'>&#10004;</font>",
        "<font color='red'>&#10006;</font>"
      )
    )
  })
  ans[["total"]] <- res
  ans
})

# *** Checkmarks for Exercise 2 ####
output$DIF_training_2_1_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_1"]])
})
output$DIF_training_2_2_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_2"]])
})
output$DIF_training_2_3_1R_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3R"]][1])
})
output$DIF_training_2_3_2R_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3R"]][2])
})
output$DIF_training_2_3_3R_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3R"]][3])
})
output$DIF_training_2_3_1F_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3F"]][1])
})
output$DIF_training_2_3_2F_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3F"]][2])
})
output$DIF_training_2_3_3F_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_3F"]][3])
})
output$DIF_training_2_4_answer <- renderUI({
  HTML(DIF_training_answers_check_2()[["check2_4"]])
})
output$DIF_training_2_answer <- renderUI({
  res <- DIF_training_answers_check_2()[["total"]]
  HTML(ifelse(is.na(res),
    "<font color = 'red'>Check the format</font>",
    ifelse(res == 1,
      "<font color='green'>Everything correct! Well done!</font>",
      paste0("<font color='red'>", round(100 * res), "% correct. Try again.</font>")
    )
  ))
})
