Likert / Top-Box Survey Analysis — Ordinal Items Done Right

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting
  #   item_names        $ named character — semantic -> humanized user names
  #   used_items        $ character — semantic item columns actually analysed
  #   dropped_df        $ data.frame(item, reason) — excluded mapped columns
  #   scale_type        $ "numeric codes" or "text labels"
  #   scale_family      $ the matched wording family, or "" for numeric
  #   common_levels     $ character — the response options, lowest first
  #   k_levels          $ integer — number of scale points
  #   top_level/bottom_level  $ the top-box and bottom-box option labels
  #   scale_consistent  $ TRUE when every item spans the same lowest/highest option
  #   scale_levels_df   $ rank, level, responses, share (pooled across items)
  #   item_summary_df   $ per item: n, top/top2/bottom box, median, mode, mean code
  #   item_dist_df      $ item, level, share (SIGNED — the diverging bar)
  #   item_rank_df      $ item, top2_box_pct (ranked, the at-a-glance bar)
  #   group_tests_df    $ per item rank-based test vs the grouping column
  #   methods_df        $ aspect, detail — full disclosure
  #   has_group / group_h / group_labels / group_test_name
  #   best_item/worst_item + their top-2-box shares
  #   mean_caveat       $ the equal-spacing sentence used wherever a mean appears
  #   n_missing_total   $ blank answers skipped across items
  #   metrics / json_output
  # === /SHARED EXPORTS ===

initial_rows <- nrow(df)
  item_cols <- grep("^item_[0-9]+$", names(df), value = TRUE)
  item_cols <- item_cols[order(as.integer(sub("^item_", "", item_cols)))]
  if (length(item_cols) < 1) {
    stop("Likert analysis needs at least one survey item mapped(item_1). Map the columns holding the scale answers.")
  }
  item_names <- setNames(humanize_semantic(item_cols, col_map), item_cols)
  all_h <- unname(item_names[item_cols])

  if (initial_rows < 10) {
    stop(sprintf("Only %d rows of survey responses were found for %s — at least 10 responses are needed to describe a response distribution.",
                 initial_rows, oxford(all_h)))
  }

dropped_items   <- character(0)
  dropped_reasons <- character(0)
  drop_item <- function(ic, reason) {
    dropped_items   <<- c(dropped_items, ic)
    dropped_reasons <<- c(dropped_reasons, reason)
  }

  raw_vals <- list(); kind <- character(0)
  for (ic in item_cols) {
    v  <- df[[ic]]
    ch <- as.character(v)
    ch[!is.na(ch) & trimws(ch) == ""] <- NA_character_
    non_blank <- !is.na(ch)
    if (sum(non_blank) == 0) {
      drop_item(ic, "every response is blank"); next
    }
    conv <- suppressWarnings(as.numeric(ch))
    if (sum(!is.na(conv)) >= 0.95 * sum(non_blank)) {
      conv[!non_blank] <- NA_real_
      raw_vals[[ic]] <- conv; kind[ic] <- "numeric"
    } else {
      raw_vals[[ic]] <- ch;   kind[ic] <- "text"
    }
  }

ok_items <- character(0); ok_rank <- list(); ok_family <- character(0)
  for (ic in names(raw_vals)) {
    v <- raw_vals[[ic]]
    obs <- v[!is.na(v)]
    lv  <- unique(obs)
    if (length(lv) < 2) {
      drop_item(ic, "every response is the same, so there is no distribution to describe"); next
    }
    if (length(lv) > 12) {
      drop_item(ic, sprintf("%d distinct answers — more than 12 distinct response options, so this reads as a measurement or an identifier rather than a rating scale",
                            length(lv)))
      next
    }
    if (kind[[ic]] == "numeric") {
      ok_items <- c(ok_items, ic)
      ok_rank[[ic]] <- sort(lv)
      ok_family[ic] <- ""
    } else {
      ord <- order_text_labels(lv)
      if (is.null(ord)) {
        drop_item(ic, "the response labels could not be placed in a meaningful order, so the answers cannot be treated as an ordered scale")
        next
      }
      ok_items <- c(ok_items, ic)
      ok_rank[[ic]] <- lv[order(ord$rank)]
      ok_family[ic] <- ord$family
    }
  }

if (length(ok_items) == 0) {
    stop(sprintf("None of the mapped columns(%s) look like a rating scale: each was blank, constant, had more than 12 distinct answers, or used labels that could not be ordered. Map the columns holding the survey scale answers.",
                 oxford(all_h)))
  }
  kinds_ok <- kind[ok_items]
  n_num <- sum(kinds_ok == "numeric"); n_txt <- sum(kinds_ok == "text")
  scale_kind <- if (n_num >= n_txt) "numeric" else "text"
  minority <- ok_items[kinds_ok != scale_kind]
  for (ic in minority) {
    drop_item(ic, sprintf("answers are %s while the other mapped items use %s, and two different answer formats cannot share one scale",
                          if (scale_kind == "numeric") "text labels" else "numeric codes",
                          if (scale_kind == "numeric") "numeric codes" else "text labels"))
  }
  used_items <- setdiff(ok_items, minority)
  scale_type <- if (scale_kind == "numeric") "numeric codes" else "text labels"

if (scale_kind == "numeric") {
    all_codes <- sort(unique(unlist(lapply(used_items, function(ic) ok_rank[[ic]]))))
    common_levels <- fmt_code(all_codes)
    level_code    <- all_codes
    scale_family  <- ""
  } else {
    lab_union <- unique(unlist(lapply(used_items, function(ic) ok_rank[[ic]])))
    ord <- order_text_labels(lab_union)
    if (is.null(ord)) {
      stop(sprintf("The mapped items(%s) use response labels from different wording families, so they cannot be placed on one shared scale. Analyse items that share a response scale together.",
                   oxford(unname(item_names[used_items]))))
    }
    common_levels <- lab_union[order(ord$rank)]
    level_code    <- rep(NA_real_, length(common_levels))
    scale_family  <- ord$family
  }
  k_levels <- length(common_levels)
  if (k_levels < 2) {
    stop(sprintf("Only one distinct response option was found across %s — a rating scale needs at least two.",
                 oxford(unname(item_names[used_items]))))
  }
  top_level    <- common_levels[k_levels]
  bottom_level <- common_levels[1]

as_index <- function(ic) {
    v <- raw_vals[[ic]]
    if (scale_kind == "numeric") match(fmt_code(v), common_levels)
    else match(as.character(v), common_levels)
  }
  idx_list <- lapply(used_items, as_index)
  names(idx_list) <- used_items

  answered_any <- Reduce(`|`, lapply(idx_list, function(z) !is.na(z)))
  final_rows   <- sum(answered_any)
  rows_removed <- initial_rows - final_rows
  if (final_rows < 10) {
    stop(sprintf("Only %d respondents answered at least one of %s — at least 10 are needed to describe a response distribution.",
                 final_rows, oxford(unname(item_names[used_items]))))
  }
  n_missing_total <- sum(vapply(idx_list, function(z) sum(is.na(z)), numeric(1)))

item_min <- vapply(idx_list, function(z) {
    zz <- z[!is.na(z)]; if (length(zz) == 0) NA_real_ else min(zz)
  }, numeric(1))
  item_max <- vapply(idx_list, function(z) {
    zz <- z[!is.na(z)]; if (length(zz) == 0) NA_real_ else max(zz)
  }, numeric(1))
  scale_consistent <- all(!is.na(item_min)) &&
    length(unique(item_min)) == 1 && length(unique(item_max)) == 1

hn <- unname(item_names[used_items])
  counts_mat <- matrix(0, nrow = length(used_items), ncol = k_levels,
                       dimnames = list(hn, common_levels))
  n_vec <- integer(length(used_items))
  median_lv <- character(length(used_items))
  mode_lv   <- character(length(used_items))
  mode_tied <- logical(length(used_items))
  mean_code <- rep(NA_real_, length(used_items))

  for (i in seq_along(used_items)) {
    z  <- idx_list[[used_items[i]]]
    zz <- z[!is.na(z)]
    n_vec[i] <- length(zz)
    cnt <- tabulate(zz, nbins = k_levels)
    counts_mat[i, ] <- cnt
    cum <- cumsum(cnt) / max(1, sum(cnt))
    hit <- which(cum >= 0.5)
    median_lv[i] <- if (length(hit) > 0) common_levels[hit[1]] else NA_character_
    best <- which(cnt == max(cnt))
    mode_lv[i]   <- common_levels[best[1]]
    mode_tied[i] <- length(best) > 1
    mean_code[i] <- if (scale_kind == "numeric") mean(level_code[zz]) else mean(zz)
  }

  share_of <- function(i, cols) 100 * sum(counts_mat[i, cols]) / max(1, n_vec[i])
  top_box  <- vapply(seq_along(used_items), function(i) share_of(i, k_levels), numeric(1))
  top2_box <- vapply(seq_along(used_items), function(i)
    share_of(i, unique(c(max(1, k_levels - 1), k_levels))), numeric(1))
  bot_box  <- vapply(seq_along(used_items), function(i) share_of(i, 1), numeric(1))

  mean_label <- if (scale_kind == "numeric") "Mean code" else "Mean position on the scale"
  mean_caveat <- paste0(
    "Reading the mean of the ", if (scale_kind == "numeric") "response codes" else "response positions",
    " as a score assumes the steps between response options are equally spaced, ",
    "which an ordinal scale does not guarantee: nothing in the data says the ",
    "distance from \"", common_levels[max(1, k_levels - 1)], "\" to \"", top_level,
    "\" matches the distance from \"", bottom_level, "\" to \"", common_levels[min(2, k_levels)],
    "\". The median category and the top-box shares use only the ordering, which is all the scale supplies."
  )

  item_summary_df <- data.frame(
    item = hn,
    responses = n_vec,
    top_box_pct = round(top_box, 1),
    top2_box_pct = round(top2_box, 1),
    bottom_box_pct = round(bot_box, 1),
    median_level = median_lv,
    mode_level = mode_lv,
    mean_code = round(mean_code, 2),
    stringsAsFactors = FALSE
  )
  ord_items <- order(-item_summary_df$top2_box_pct, item_summary_df$item)
  item_summary_df <- item_summary_df[ord_items, , drop = FALSE]
  rownames(item_summary_df) <- NULL

  best_item  <- item_summary_df$item[1]
  best_top2  <- item_summary_df$top2_box_pct[1]
  worst_item <- item_summary_df$item[nrow(item_summary_df)]
  worst_top2 <- item_summary_df$top2_box_pct[nrow(item_summary_df)]

  item_rank_df <- data.frame(
    item = item_summary_df$item,
    top2_box_pct = item_summary_df$top2_box_pct,
    stringsAsFactors = FALSE
  )

pooled <- colSums(counts_mat)
  scale_levels_df <- data.frame(
    rank = seq_len(k_levels),
    level = common_levels,
    responses = as.integer(pooled),
    share = round(100 * pooled / max(1, sum(pooled)), 1),
    stringsAsFactors = FALSE
  )

midpoint <- (k_levels + 1) / 2
  sign_of  <- ifelse(seq_len(k_levels) <= midpoint, -1, 1)
  mid_on_left <- (k_levels %% 2) == 1
  dist_rows <- list()
  for (i in seq_along(used_items)) {
    ii <- match(hn[i], item_summary_df$item)
    for (j in seq_len(k_levels)) {
      share <- 100 * counts_mat[i, j] / max(1, n_vec[i])
      dist_rows[[length(dist_rows) + 1]] <- data.frame(
        item = hn[i], level = common_levels[j],
        share = round(sign_of[j] * share, 2),
        order_key = ii, level_rank = j,
        stringsAsFactors = FALSE
      )
    }
  }
  item_dist_df <- do.call(rbind, dist_rows)
  item_dist_df <- item_dist_df[order(item_dist_df$order_key, item_dist_df$level_rank), , drop = FALSE]
  item_dist_df <- item_dist_df[, c("item", "level", "share")]
  rownames(item_dist_df) <- NULL

group_h <- humanize_semantic("group", col_map)
  has_group <- "group" %in% names(df)
  group_labels <- character(0); group_test_name <- ""
  group_note <- ""
  gvec <- NULL
  if (has_group) {
    g <- as.character(df$group)
    g[!is.na(g) & trimws(g) == ""] <- NA_character_
    tab <- table(g)
    small <- names(tab)[tab < 5]
    if (length(small) > 0) g[g %in% small] <- NA_character_
    tab <- sort(table(g), decreasing = TRUE)
    if (length(tab) > 8) {
      keep <- names(tab)[1:8]
      g[!is.na(g) & !(g %in% keep)] <- NA_character_
      group_note <- sprintf("Only the 8 largest %s values were compared; smaller ones were set aside. ", group_h)
    }
    if (length(small) > 0) {
      group_note <- paste0(group_note, sprintf(
        "%d %s value(s) with fewer than 5 responses were set aside as too small to compare. ",
        length(small), group_h))
    }
    gvec <- g
    group_labels <- sort(unique(g[!is.na(g)]))
    if (length(group_labels) < 2) {
      has_group <- FALSE
      group_note <- paste0(group_note, sprintf(
        "Fewer than two %s values had enough responses, so no comparison was run. ", group_h))
    } else {
      group_test_name <- if (length(group_labels) == 2)
        "Mann-Whitney U(Wilcoxon rank-sum)" else "Kruskal-Wallis"
    }
  }

  if (has_group) {
    rows <- list(); praw <- numeric(0)
    for (i in seq_along(used_items)) {
      z <- idx_list[[used_items[i]]]
      ok <- !is.na(z) & !is.na(gvec)
      zz <- z[ok]; gg <- factor(gvec[ok])
      stat <- NA_real_; pv <- NA_real_; eff <- NA_real_
      hi <- "n/a"; lo <- "n/a"; tname <- group_test_name
      if (length(zz) >= 10 && nlevels(gg) >= 2 &&
          all(table(gg) >= 3) && length(unique(zz)) > 1) {
        mr <- tapply(zz, gg, mean)
        mr_ok <- which(!is.na(mr))
        if (length(mr_ok) > 0) {
          hi <- names(mr)[mr_ok[which.max(mr[mr_ok])]]
          lo <- names(mr)[mr_ok[which.min(mr[mr_ok])]]
        }
        if (nlevels(gg) == 2) {
          tt <- suppressWarnings(tryCatch(stats::wilcox.test(zz ~ gg), error = function(e) NULL))
          if (!is.null(tt)) {
            n1 <- sum(gg == levels(gg)[1]); n2 <- sum(gg == levels(gg)[2])
            stat <- unname(tt$statistic); pv <- tt$p.value
            eff <- 2 * stat / (n1 * n2) - 1
          }
        } else {
          tt <- tryCatch(stats::kruskal.test(zz ~ gg), error = function(e) NULL)
          if (!is.null(tt)) {
            stat <- unname(tt$statistic); pv <- tt$p.value
            eff <- stat / max(1, length(zz) - 1)
          }
        }
      } else {
        tname <- "not run"
      }
      praw <- c(praw, pv)
      rows[[i]] <- list(item = hn[i], test = tname, statistic = stat,
                        p_value = pv, effect = eff, hi = hi, lo = lo,
                        n = length(zz))
    }
    padj <- rep(NA_real_, length(praw))
    okp <- which(!is.na(praw))
    if (length(okp) > 0) padj[okp] <- stats::p.adjust(praw[okp], method = "holm")

    eff_label <- if (length(group_labels) == 2) "rank-biserial correlation" else "epsilon squared"
    group_tests_df <- do.call(rbind, lapply(seq_along(rows), function(i) {
      r <- rows[[i]]
      sig <- !is.na(padj[i]) && padj[i] < 0.05
      interp <- if (r$test == "not run") {
        sprintf("Too few usable responses to compare %s across %s.", r$item, group_h)
      } else if (sig) {
        sprintf("%s answer %s differently across %s(%s after Holm correction across %d items): %s rates it highest and %s lowest. The test compares the ORDER of the answers, so it needs no assumption about the spacing between response options.",
                r$item, "is answered", group_h, fmt_pp(padj[i]), length(used_items), r$hi, r$lo)
      } else {
        sprintf("No difference in how %s is answered across %s survives correction(%s after Holm correction across %d items).",
                r$item, group_h, fmt_pp(padj[i]), length(used_items))
      }
      data.frame(item = r$item, test = r$test,
                 statistic = if (is.na(r$statistic)) NA_real_ else round(r$statistic, 3),
                 p_value = fmt_p(r$p_value), p_adjusted = fmt_p(padj[i]),
                 effect_size = if (is.na(r$effect)) NA_real_ else round(r$effect, 3),
                 highest_group = r$hi, lowest_group = r$lo,
                 interpretation = interp, stringsAsFactors = FALSE)
    }))
    rownames(group_tests_df) <- NULL
    group_sig_items <- group_tests_df$item[!is.na(padj) & padj < 0.05]
    n_group_sig <- length(group_sig_items)
  } else {
    eff_label <- "n/a"
    group_tests_df <- data.frame(
      item = "(no comparison run)", test = "not run", statistic = NA_real_,
      p_value = "n/a", p_adjusted = "n/a", effect_size = NA_real_,
      highest_group = "n/a", lowest_group = "n/a",
      interpretation = paste0(
        if (nzchar(group_note)) group_note else
          sprintf("No %s column was mapped, so the responses were described for the whole sample only. ", group_h),
        "Map a categorical column to compare response distributions between groups with a rank-based test."),
      stringsAsFactors = FALSE)
    group_sig_items <- character(0); n_group_sig <- 0
  }

dropped_df <- if (length(dropped_items) > 0) {
    data.frame(item = unname(item_names[dropped_items]),
               reason = dropped_reasons, stringsAsFactors = FALSE)
  } else {
    data.frame(item = character(0), reason = character(0), stringsAsFactors = FALSE)
  }

  scale_sentence <- if (scale_kind == "numeric") {
    sprintf("A %d-point scale was detected from the data: the response codes %s, read in ascending order.",
            k_levels, paste(common_levels, collapse = ", "))
  } else {
    sprintf("A %d-point scale was detected from the data: the labels %s, ordered as a %s scale from the wording itself rather than assumed.",
            k_levels, paste0("\"", paste(common_levels, collapse = "\", \""), "\""), scale_family)
  }

  methods_df <- data.frame(
    aspect = c("Detected scale", "Top box", "Top-2 box", "Bottom box",
               "Median category", "Modal category", mean_label,
               "Diverging bar", "Group comparison", "Missing answers"),
    detail = c(
      scale_sentence,
      sprintf("The share of answers at the highest option, \"%s\".", top_level),
      sprintf("The share of answers at the two highest options, \"%s\" and \"%s\".",
              common_levels[max(1, k_levels - 1)], top_level),
      sprintf("The share of answers at the lowest option, \"%s\".", bottom_level),
      "The lowest option whose cumulative share of answers reaches 50 percent. It is defined for ties and for an even number of respondents alike, and it uses only the ordering of the options.",
      "The single most-chosen option. Where two options tie, the lower of the two is reported and the tie is flagged.",
      mean_caveat,
      sprintf("Shares at or below the scale midpoint are drawn left of zero and shares above it to the right, with items ordered by their top-2-box share%s.",
              if (mid_on_left) sprintf("; the middle option \"%s\" is drawn on the left so the right-hand side reads as the top half of the scale only",
                                       common_levels[ceiling(midpoint)]) else ""),
      if (has_group)
        sprintf("%s on the ranked answers, one test per item, %s correction across the %d items, with %s as the effect size. A rank-based test is used rather than a t-test on the raw codes because the codes are labels in order, not measured quantities.",
                group_test_name, "Holm", length(used_items), eff_label)
      else
        group_tests_df$interpretation[1],
      sprintf("%s blank answers were skipped; each item&#x27;s percentages use that item's own answered count as the denominator, so items with different response rates stay comparable.",
              comma(n_missing_total))
    ),
    stringsAsFactors = FALSE
  )

  metrics <- list(
    `Items Analysed`      = length(used_items),
    `Respondents`         = final_rows,
    `Scale Points`        = k_levels,
    `Scale Type`          = scale_type,
    `Top Box`             = top_level,
    `Highest Top-2-Box`   = paste0(best_item, " (", pct1(best_top2), ")"),
    `Lowest Top-2-Box`    = paste0(worst_item, " (", pct1(worst_top2), ")"),
    `Group Differences`   = if (!has_group) "not tested" else
      paste0(n_group_sig, " of ", length(used_items), " items")
  )

  group_clause <- if (has_group) {
    if (n_group_sig > 0)
      paste0(" Across ", group_h, ", ", oxford(group_sig_items),
             " differ(s) by a ", group_test_name,
             " test after Holm correction; the remaining ",
             length(used_items) - n_group_sig, " do not.")
    else
      paste0(" No item differs across ", group_h,
             " once the ", group_test_name,
             " results are Holm-corrected across the ", length(used_items), " items.")
  } else ""

  json_output <- list(
    answer = paste0(
      "Ordinal analysis of ", length(used_items), " survey item(s) answered by ",
      comma(final_rows), " respondents on a ", k_levels, "-point scale of ",
      scale_type, " (", paste(common_levels, collapse = " < "), "). ",
      best_item, " scores highest with ", pct1(best_top2),
      " in the top two options and ", worst_item, " lowest with ",
      pct1(worst_top2), "; the top box is \"", top_level, "\".",
      group_clause,
      " Percentages and median categories use only the ordering of the response options; ",
      "where a mean of the codes is reported it assumes equal spacing between options, which the scale does not guarantee."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "scale_detection", "item_summary",
        "diverging_distribution", "top_box_ranking", "group_comparison", "methods"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows, rows_removed = rows_removed,
    item_names = item_names, used_items = used_items, item_h = hn,
    dropped_df = dropped_df, n_missing_total = n_missing_total,
    scale_type = scale_type, scale_kind = scale_kind, scale_family = scale_family,
    common_levels = common_levels, k_levels = k_levels,
    top_level = top_level, bottom_level = bottom_level,
    scale_consistent = scale_consistent, scale_sentence = scale_sentence,
    mid_on_left = mid_on_left, mean_caveat = mean_caveat, mean_label = mean_label,
    scale_levels_df = scale_levels_df, item_summary_df = item_summary_df,
    item_dist_df = item_dist_df, item_rank_df = item_rank_df,
    group_tests_df = group_tests_df, methods_df = methods_df,
    has_group = has_group, group_h = group_h, group_labels = group_labels,
    group_test_name = group_test_name, group_note = group_note,
    group_sig_items = group_sig_items, n_group_sig = n_group_sig,
    mode_tied_any = any(mode_tied),
    best_item = best_item, best_top2 = best_top2,
    worst_item = worst_item, worst_top2 = worst_top2,
    metrics = metrics, json_output = json_output
  )
}