Conjoint Analysis — Feature and Price Trade-offs

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
  #   shape              $ "choice" or "ratings" — which data shape was detected
  #   shape_reason       $ computed sentence explaining the detection
  #   pref_h/task_h/resp_h $ humanized user names for the mapped singles
  #   attr_h             $ named character — semantic attribute -> user name
  #   used_attrs         $ character — semantic attribute columns actually used
  #   dropped_df         $ data.frame(column, reason) — excluded attributes
  #   levels_list        $ named list — the tested levels per used attribute
  #   pw_df              $ part-worth table (attribute, level, utility, CI, ...)
  #   imp_df             $ attribute importance table incl. the tested range
  #   wtp_df             $ willingness-to-pay table (or a single reason row)
  #   wtp_ok             $ TRUE when WTP was computed
  #   sim_df             $ market-simulation table
  #   het_df             $ heterogeneity table (or a single reason row)
  #   het_ok / het_flag  $ whether it ran / whether a split audience is implied
  #   methods_df         $ item/detail disclosure table
  #   fit_stat/fit_label $ the fit statistic and its name
  #   n_resp/n_tasks     $ respondents and choice tasks (choice shape)
  #   price_attr/price_slope/price_slope_ci  $ price detection results
  #   stated_note        $ the stated-preference caveat, used in several cards
  #   metrics / json_output
  # === /SHARED EXPORTS ===

  initial_rows <- nrow(df)

pref_h <- humanize_semantic("preference", col_map)
  task_h <- humanize_semantic("task", col_map)
  resp_h <- humanize_semantic("respondent", col_map)

  attr_cols <- grep("^attribute_[0-9]+$", names(df), value = TRUE)
  attr_cols <- attr_cols[order(as.integer(sub("^attribute_", "", attr_cols)))]
  attr_h <- setNames(humanize_semantic(attr_cols, col_map), attr_cols)

  if (!("preference" %in% names(df))) {
    stop(sprintf("Conjoint analysis needs the preference column &#x27;%s' mapped — the rating, rank score, or 0/1 chosen flag that each row carries.", pref_h))
  }
  if (length(attr_cols) < 1) {
    stop("Conjoint analysis needs at least one product attribute column mapped(attribute_1, attribute_2, ...) — the columns holding each profile&#x27;s feature levels.")
  }

  has_task <- "task" %in% names(df)
  has_resp <- "respondent" %in% names(df)

pv_raw <- df$preference
  if (is.logical(pv_raw)) {
    pref <- as.numeric(pv_raw)
  } else if (is.numeric(pv_raw)) {
    pref <- as.numeric(pv_raw)
  } else {
    ch <- as.character(pv_raw)
    non_blank <- !is.na(ch) & trimws(ch) != ""
    conv <- suppressWarnings(as.numeric(ch))
    if (sum(non_blank) == 0 ||
        sum(!is.na(conv[non_blank])) < 0.95 * sum(non_blank)) {
      stop(sprintf("The preference column &#x27;%s' does not look numeric — fewer than 95%% of its values parse as numbers. Map a rating, a rank score, or a 0/1 chosen flag.", pref_h))
    }
    pref <- conv
  }

drop_col <- character(0); drop_reason <- character(0)
  attrs <- list(); levels_list <- list(); lumped <- character(0)
  MAX_LEVELS_KEPT <- 12
  for (ac in attr_cols) {
    lv_raw <- trimws(as.character(df[[ac]]))
    lv_raw[is.na(lv_raw) | lv_raw == ""] <- "Missing"
    nd <- length(unique(lv_raw))
    if (nd < 2) {
      drop_col <- c(drop_col, attr_h[[ac]])
      drop_reason <- c(drop_reason, sprintf("constant — every row has the same value(&#x27;%s'), so it carries no trade-off information", unique(lv_raw)[1]))
      next
    }
    if (nd >= 0.5 * initial_rows || nd > 20) {
      drop_col <- c(drop_col, attr_h[[ac]])
      drop_reason <- c(drop_reason, sprintf("%s distinct values across %s rows — that is an identifier or free text, not a designed attribute with a small set of tested levels", ct(nd), ct(initial_rows)))
      next
    }
    if (nd > MAX_LEVELS_KEPT) {
      tab <- sort(table(lv_raw), decreasing = TRUE)
      keep <- names(tab)[seq_len(MAX_LEVELS_KEPT - 1)]
      lv_raw[!(lv_raw %in% keep)] <- "Other"
      lumped <- c(lumped, attr_h[[ac]])
    }
    attrs[[ac]] <- lv_raw
    levels_list[[ac]] <- sort(unique(lv_raw))
  }
  used_attrs <- names(attrs)
  if (length(used_attrs) < 1) {
    stop(sprintf("None of the mapped attribute columns(%s) can be used: each is either constant or looks like an identifier rather than a designed attribute with a small set of tested levels.",
                 paste(unname(attr_h), collapse = ", ")))
  }
  dropped_df <- if (length(drop_col) > 0) {
    data.frame(column = drop_col, reason = drop_reason, stringsAsFactors = FALSE)
  } else {
    data.frame(column = character(0), reason = character(0), stringsAsFactors = FALSE)
  }

A <- as.data.frame(attrs, stringsAsFactors = FALSE)
  names(A) <- used_attrs
  keep <- !is.na(pref)
  for (ac in used_attrs) keep <- keep & !is.na(A[[ac]])
  pref <- pref[keep]
  A <- A[keep, , drop = FALSE]
  task_vec <- if (has_task) as.character(df$task)[keep] else NULL
  resp_vec <- if (has_resp) as.character(df$respondent)[keep] else NULL
  n_incomplete <- sum(!keep)

for (ac in used_attrs) levels_list[[ac]] <- sort(unique(A[[ac]]))
  still_varying <- vapply(used_attrs, function(ac) length(levels_list[[ac]]) >= 2, TRUE)
  if (any(!still_varying)) {
    lost <- used_attrs[!still_varying]
    dropped_df <- rbind(dropped_df, data.frame(
      column = unname(attr_h[lost]),
      reason = "became constant once rows with a missing preference value were dropped",
      stringsAsFactors = FALSE))
    used_attrs <- used_attrs[still_varying]
    A <- A[, used_attrs, drop = FALSE]
    levels_list <- levels_list[used_attrs]
  }
  if (length(used_attrs) < 1) {
    stop(sprintf("No attribute column still varies once rows with a missing &#x27;%s' value are dropped, so no part-worth utilities can be estimated.", pref_h))
  }

n_kept <- length(pref)
  uniq_pref <- sort(unique(pref))
  binary_pref <- length(uniq_pref) == 2 && all(uniq_pref %in% c(0, 1))
  shape <- "ratings"; shape_reason <- ""; n_tasks <- NA_integer_
  bad_tasks <- 0L
  if (has_task && binary_pref) {
    tsize <- table(task_vec)
    tchosen <- tapply(pref, task_vec, sum)
    good <- names(tsize)[tsize >= 2 & !is.na(tchosen[names(tsize)]) &
                           tchosen[names(tsize)] == 1]
    bad_tasks <- as.integer(length(tsize) - length(good))
    if (length(good) >= 20) {
      sel <- task_vec %in% good
      pref <- pref[sel]; A <- A[sel, , drop = FALSE]
      task_vec <- task_vec[sel]
      if (!is.null(resp_vec)) resp_vec <- resp_vec[sel]
      for (ac in used_attrs) levels_list[[ac]] <- sort(unique(A[[ac]]))
      shape <- "choice"
      n_tasks <- length(good)
      shape_reason <- sprintf(
        "&#x27;%s' takes only the values 0 and 1 and '%s' groups the rows into choice sets, so the data was read as choice-based conjoint: %s choice tasks, each with the alternatives that were shown together and exactly one chosen.",
        pref_h, task_h, ct(n_tasks))
    } else {
      shape_reason <- sprintf(
        "&#x27;%s' is a 0/1 flag and '%s' was mapped, but only %s choice tasks have at least two alternatives and exactly one chosen — too few for a conditional logit, so the data was read as ratings-based instead.",
        pref_h, task_h, ct(length(good)))
    }
  }
  if (shape == "ratings" && shape_reason == "") {
    shape_reason <- if (binary_pref && !has_task) {
      sprintf("&#x27;%s' takes only the values 0 and 1 but no choice-task column was mapped, so the alternatives that competed with each other are unknown. The data was fitted as a linear probability model on the 0/1 flag; map the choice-task column to get a proper choice-based conjoint.", pref_h)
    } else {
      sprintf("&#x27;%s' carries %s distinct values, so the data was read as a ratings or rankings conjoint: one row per profile evaluated, fitted by ordinary least squares on dummy-coded attribute levels.",
              pref_h, ct(length(uniq_pref)))
    }
  }

cols <- list(); meta_attr <- character(0); meta_level <- character(0)
  for (ac in used_attrs) {
    lv <- levels_list[[ac]]
    for (k in seq_along(lv)[-1]) {
      cols[[length(cols) + 1L]] <- as.numeric(A[[ac]] == lv[k])
      meta_attr <- c(meta_attr, ac); meta_level <- c(meta_level, lv[k])
    }
  }
  X <- do.call(cbind, cols)
  colnames(X) <- paste0("p", seq_len(ncol(X)))
  n_params <- ncol(X)
  final_rows <- nrow(A)
  rows_removed <- initial_rows - final_rows
  min_rows <- max(20L, as.integer(3 * n_params))
  if (final_rows < min_rows) {
    stop(sprintf("Only %s usable rows of &#x27;%s' across the mapped attributes (%s) remained — a conjoint fit for %s attribute levels needs at least %s rows.",
                 ct(final_rows), pref_h,
                 paste(unname(attr_h[used_attrs]), collapse = ", "),
                 ct(n_params + length(used_attrs)), ct(min_rows)))
  }

fit_label <- ""; fit_stat <- NA_real_; resid_sd <- NA_real_
  base_utility <- NA_real_; converged <- TRUE; ll_final <- NA_real_
  if (shape == "ratings") {
    fdat <- as.data.frame(X); fdat$y__ <- pref
    fit <- stats::lm(y__ ~ ., data = fdat)
    b <- stats::coef(fit)
    if (any(is.na(b))) {
      alias <- names(b)[is.na(b)]
      bad <- unique(unname(attr_h[meta_attr[match(alias, colnames(X))]]))
      stop(sprintf("Some attribute levels cannot be separated from the others in this design(%s) — the profiles shown never varied them independently, so their part-worths are not identified.",
                   paste(bad[!is.na(bad)], collapse = ", ")))
    }
    Vfull <- stats::vcov(fit)
    b_slope <- b[-1]; V <- Vfull[-1, -1, drop = FALSE]
    fit_stat <- summary(fit)$r.squared
    fit_label <- "R-squared"
    resid_sd <- stats::sd(stats::residuals(fit))
    intercept <- unname(b[1])
  } else {
    ord <- order(task_vec)
    X <- X[ord, , drop = FALSE]
    pref <- pref[ord]; y <- pref
    tv <- task_vec[ord]
    if (!is.null(resp_vec)) resp_vec <- resp_vec[ord]
    A <- A[ord, , drop = FALSE]
    tid <- as.integer(factor(tv, levels = unique(tv)))
    negll <- function(bb) {
      eta <- as.vector(X %*% bb)
      if (any(!is.finite(eta))) return(1e10)
      eta <- eta - mean(eta)
      if (max(abs(eta)) > 400) return(1e10)
      e <- exp(eta)
      d <- as.vector(rowsum(e, tid, reorder = FALSE))[tid]
      -sum(y * (eta - log(d)))
    }
    gradll <- function(bb) {
      eta <- as.vector(X %*% bb)
      if (any(!is.finite(eta))) return(rep(0, length(bb)))
      eta <- eta - mean(eta)
      if (max(abs(eta)) > 400) return(rep(0, length(bb)))
      e <- exp(eta)
      d <- as.vector(rowsum(e, tid, reorder = FALSE))[tid]
      -as.vector(crossprod(X, y - e / d))
    }
    opt <- stats::optim(rep(0, n_params), negll, gradll, method = "BFGS",
                        control = list(maxit = 500, reltol = 1e-12))
    converged <- isTRUE(opt$convergence == 0)
    b_slope <- opt$par; names(b_slope) <- colnames(X)
    ll_final <- -opt$value
    eta <- as.vector(X %*% b_slope); eta <- eta - mean(eta)
    e <- exp(eta)
    d <- as.vector(rowsum(e, tid, reorder = FALSE))[tid]
    p_hat <- e / d

W <- X * p_hat
    Info <- crossprod(X, W) - crossprod(rowsum(W, tid, reorder = FALSE))
    V <- tryCatch(solve(Info), error = function(e) NULL)
    if (is.null(V)) {
      stop(sprintf("The choice design does not identify all attribute levels(%s) — some levels never varied independently within a choice task, so their part-worths cannot be separated.",
                   paste(unname(attr_h[used_attrs]), collapse = ", ")))
    }

tsz <- as.vector(rowsum(rep(1, length(y)), tid, reorder = FALSE))
    ll_null <- -sum(log(tsz))
    fit_stat <- 1 - ll_final / ll_null
    fit_label <- "McFadden Pseudo R-squared"
    intercept <- NA_real_
    final_rows <- nrow(A)
    rows_removed <- initial_rows - final_rows
  }

lvl_attr <- character(0); lvl_name <- character(0)
  Mrows <- list()
  for (ac in used_attrs) {
    lv <- levels_list[[ac]]
    K <- length(lv)
    idx <- which(meta_attr == ac)          # parameter positions, levels 2..K
    for (k in seq_len(K)) {
      row <- rep(0, n_params)
      row[idx] <- -1 / K
      if (k > 1) row[idx[k - 1]] <- row[idx[k - 1]] + 1
      Mrows[[length(Mrows) + 1L]] <- row
      lvl_attr <- c(lvl_attr, ac); lvl_name <- c(lvl_name, lv[k])
    }
  }
  M <- do.call(rbind, Mrows)
  u <- as.vector(M %*% b_slope)
  Vu <- M %*% V %*% t(M)
  u_se <- sqrt(pmax(0, diag(Vu)))
  zq <- 1.959964
  pw_df <- data.frame(
    attribute = unname(attr_h[lvl_attr]),
    level = lvl_name,
    level_label = paste0(unname(attr_h[lvl_attr]), ": ", lvl_name),
    utility = round(u, 4),
    ci_low = round(u - zq * u_se, 4),
    ci_high = round(u + zq * u_se, 4),
    std_error = round(u_se, 4),
    stringsAsFactors = FALSE
  )
  if (shape == "ratings") {

shift <- 0
    for (ac in used_attrs) {
      idx <- which(meta_attr == ac); K <- length(levels_list[[ac]])
      shift <- shift + sum(b_slope[idx]) / K
    }
    base_utility <- intercept + shift
  }

imp_attr <- character(0); imp_range <- numeric(0)
  imp_best <- character(0); imp_worst <- character(0); imp_tested <- character(0)
  for (ac in used_attrs) {
    sel <- lvl_attr == ac
    uu <- u[sel]; nm <- lvl_name[sel]
    hi <- safe_which_max(uu); lo <- safe_which_min(uu)
    imp_attr <- c(imp_attr, unname(attr_h[[ac]]))
    imp_range <- c(imp_range, if (is.na(hi) || is.na(lo)) NA_real_ else uu[hi] - uu[lo])
    imp_best <- c(imp_best, if (is.na(hi)) "n/a" else nm[hi])
    imp_worst <- c(imp_worst, if (is.na(lo)) "n/a" else nm[lo])
    imp_tested <- c(imp_tested, paste(nm, collapse = " | "))
  }
  tot_range <- sum(imp_range, na.rm = TRUE)
  imp_pct <- if (tot_range > 0) 100 * imp_range / tot_range else rep(NA_real_, length(imp_range))
  ord_i <- order(-imp_pct, na.last = TRUE)
  imp_df <- data.frame(
    attribute = imp_attr[ord_i],
    importance_pct = round(imp_pct[ord_i], 2),
    utility_range = round(imp_range[ord_i], 4),
    best_level = imp_best[ord_i],
    worst_level = imp_worst[ord_i],
    tested_levels = imp_tested[ord_i],
    n_levels = as.integer(vapply(used_attrs[ord_i], function(ac) length(levels_list[[ac]]), 1L)),
    stringsAsFactors = FALSE
  )
  top_imp <- imp_df[1, ]

price_attr <- NA_character_; price_sym <- ""
  price_levels_num <- NULL; price_idx <- NULL
  for (ac in used_attrs) {
    lbl <- levels_list[[ac]]
    hname <- unname(attr_h[[ac]])
    name_hit <- grepl("pric|cost|fee|charg|tariff|msrp|rate|\\$|usd|eur|gbp", hname, ignore.case = TRUE)
    sym_hit <- all(grepl("^[\\$€£¥]", lbl))
    if (!name_hit && !sym_hit) next
    nums <- suppressWarnings(as.numeric(gsub("[^0-9.-]", "", lbl)))
    if (any(is.na(nums)) || length(unique(nums)) < 2) next
    price_attr <- ac
    price_idx <- which(lvl_attr == ac)
    price_levels_num <- nums[match(lvl_name[price_idx], lbl)]
    syms <- substr(lbl, 1, 1)
    price_sym <- if (length(unique(syms)) == 1 && grepl("[\\$€£¥]", syms[1])) syms[1] else ""
    break
  }
  price_slope <- NA_real_; price_slope_se <- NA_real_
  price_slope_lo <- NA_real_; price_slope_hi <- NA_real_
  price_lin_r2 <- NA_real_; wtp_ok <- FALSE; wtp_reason <- ""
  cvec <- NULL
  if (!is.na(price_attr)) {
    pnum <- price_levels_num
    dev <- pnum - mean(pnum); ssd <- sum(dev^2)
    cvec <- rep(0, length(u)); cvec[price_idx] <- dev / ssd
    price_slope <- sum(cvec * u)
    price_slope_se <- sqrt(max(0, as.numeric(t(cvec) %*% Vu %*% cvec)))
    price_slope_lo <- price_slope - zq * price_slope_se
    price_slope_hi <- price_slope + zq * price_slope_se
    up <- u[price_idx]
    fitted_up <- mean(up) + price_slope * dev
    sst <- sum((up - mean(up))^2)
    price_lin_r2 <- if (sst > 0) 1 - sum((up - fitted_up)^2) / sst else NA_real_
    if (is.na(price_slope) || price_slope >= 0) {
      wtp_reason <- sprintf("The fitted price coefficient is %s per unit of &#x27;%s' — not negative, so respondents did not trade utility away as price rose in this data. Dividing by it would produce a nonsense willingness-to-pay, so none is reported.",
                            r3(price_slope), unname(attr_h[[price_attr]]))
    } else if (price_slope_hi >= 0) {
      wtp_reason <- sprintf("The price coefficient is %s per unit of &#x27;%s', but its 95%% interval (%s to %s) still includes zero, so the exchange rate between utility and money is not established. Willingness-to-pay is a division by that coefficient and would be unbounded here, so none is reported.",
                            r3(price_slope), unname(attr_h[[price_attr]]),
                            r3(price_slope_lo), r3(price_slope_hi))
    } else {
      wtp_ok <- TRUE
    }
  } else {
    wtp_reason <- sprintf("No price attribute could be identified among the mapped attributes(%s): willingness-to-pay needs one attribute whose levels are money amounts. Map the price column as one of the attributes to get it.",
                          paste(unname(attr_h[used_attrs]), collapse = ", "))
  }

money_per_util <- NA_real_
  if (wtp_ok) {
    dslope <- -price_slope                     # positive utility cost per unit
    money_per_util <- 1 / dslope
    cov_u_s <- as.vector(Vu %*% cvec)          # Cov(u_l, slope) for every level
    keepw <- setdiff(seq_along(u), price_idx)
    wtp <- u[keepw] / dslope
    var_wtp <- (1 / dslope^2) * diag(Vu)[keepw] +
      (u[keepw]^2 / dslope^4) * price_slope_se^2 +
      (2 * u[keepw] / dslope^3) * cov_u_s[keepw]
    wtp_se <- sqrt(pmax(0, var_wtp))
    o <- order(-wtp)
    wtp_df <- data.frame(
      level_label = pw_df$level_label[keepw][o],
      wtp = round(wtp[o], 2),
      wtp_low = round((wtp - zq * wtp_se)[o], 2),
      wtp_high = round((wtp + zq * wtp_se)[o], 2),
      utility = round(u[keepw][o], 4),
      stringsAsFactors = FALSE
    )
    wtp_df$interpretation <- sprintf(
      "Stated-preference estimate: relative to the average level of its own attribute, this level is worth about %s%s in the units of &#x27;%s'. It is not a price a customer has agreed to pay.",
      price_sym, r2(abs(wtp_df$wtp)), unname(attr_h[[price_attr]]))
    wtp_df$interpretation[wtp_df$wtp < 0] <- sprintf(
      "Stated-preference estimate: relative to the average level of its own attribute, this level costs about %s%s of value in the units of &#x27;%s'. It is not a price a customer has agreed to pay.",
      price_sym, r2(abs(wtp_df$wtp[wtp_df$wtp < 0])), unname(attr_h[[price_attr]]))
  } else {
    wtp_df <- data.frame(
      level_label = "Willingness-to-pay not reported",
      wtp = NA_real_, wtp_low = NA_real_, wtp_high = NA_real_,
      utility = NA_real_, interpretation = wtp_reason,
      stringsAsFactors = FALSE
    )
  }

best_lv <- character(0); worst_lv <- character(0); modal_lv <- character(0)
  for (ac in used_attrs) {
    sel <- lvl_attr == ac
    uu <- u[sel]; nm <- lvl_name[sel]
    hi <- safe_which_max(uu); lo <- safe_which_min(uu)
    best_lv <- c(best_lv, if (is.na(hi)) nm[1] else nm[hi])
    worst_lv <- c(worst_lv, if (is.na(lo)) nm[1] else nm[lo])
    tb <- table(A[[ac]])
    modal_lv <- c(modal_lv, names(tb)[safe_which_max(as.numeric(tb))])
  }
  names(best_lv) <- names(worst_lv) <- names(modal_lv) <- used_attrs
  util_of <- function(cfg) {
    s <- 0
    for (ac in used_attrs) {
      j <- which(lvl_attr == ac & lvl_name == cfg[[ac]])
      if (length(j) == 1) s <- s + u[j]
    }
    s
  }
  cfg_list <- list(); cfg_names <- character(0)
  cfg_list[[1]] <- best_lv;  cfg_names[1] <- "Best level on every attribute"
  cfg_list[[2]] <- modal_lv; cfg_names[2] <- "Most common configuration in the data"
  cfg_list[[3]] <- worst_lv; cfg_names[3] <- "Worst level on every attribute"
  if (!is.na(price_attr) && length(price_levels_num) >= 2) {
    pnm <- lvl_name[price_idx]
    hi_price <- pnm[safe_which_max(price_levels_num)]
    lo_price <- pnm[safe_which_min(price_levels_num)]
    c4 <- best_lv;  c4[[price_attr]] <- hi_price
    c5 <- modal_lv; c5[[price_attr]] <- lo_price
    cfg_list[[4]] <- c4
    cfg_names[4] <- sprintf("Best features at the highest tested %s(%s)",
                            unname(attr_h[[price_attr]]), hi_price)
    cfg_list[[5]] <- c5
    cfg_names[5] <- sprintf("Most common features at the lowest tested %s(%s)",
                            unname(attr_h[[price_attr]]), lo_price)
  }
  spec_str <- vapply(cfg_list, function(cfg)
    paste(sprintf("%s = %s", unname(attr_h[used_attrs]),
                  unlist(cfg[used_attrs])), collapse = "; "), "")
  dup <- duplicated(spec_str)
  cfg_list <- cfg_list[!dup]; cfg_names <- cfg_names[!dup]; spec_str <- spec_str[!dup]
  cfg_util <- vapply(cfg_list, util_of, 0)
  if (shape == "choice") {
    ee <- exp(cfg_util - max(cfg_util))
    cfg_share <- round(100 * ee / sum(ee), 2)
    share_basis <- "Logit choice share among the configurations listed here, on the scale the choice model identified."
  } else {
    cfg_share <- rep(NA_real_, length(cfg_util))
    share_basis <- sprintf("Not identified: turning aggregate &#x27;%s' utilities into choice shares needs an arbitrary scale factor that ratings data does not pin down, so the ranking and the utility gaps are reported instead of invented shares.", pref_h)
  }
  ord_s <- order(-cfg_util)
  sim_df <- data.frame(
    configuration = cfg_names[ord_s],
    total_utility = round(cfg_util[ord_s], 4),
    predicted_share_pct = cfg_share[ord_s],
    levels_used = spec_str[ord_s],
    share_basis = share_basis,
    stringsAsFactors = FALSE
  )

het_ok <- FALSE; het_flag <- FALSE; het_reason <- ""
  het_df <- NULL; het_min <- NA_real_; het_worst_attr <- NA_character_
  n_resp <- if (!is.null(resp_vec)) length(unique(resp_vec)) else NA_integer_
  if (shape == "choice") {
    het_reason <- sprintf("The respondent-level check is not run on choice data. Each respondent contributes only a handful of binary choices, so their individual part-worths are not estimable without a hierarchical model, which this module deliberately does not fit rather than fit badly. Read every aggregate part-worth near zero as ambiguous: it may mean indifference, or it may mean one group of respondents wants the level and another rejects it.")
  } else if (is.null(resp_vec)) {
    het_reason <- "No respondent column was mapped, so the aggregate part-worths cannot be checked for a split audience. Map the respondent identifier to find out whether a level averaging to neutral is genuinely neutral or is loved by one group and disliked by another."
  } else if (n_resp < 8) {
    het_reason <- sprintf("Only %s respondents were identified in &#x27;%s' — too few to check whether the aggregate part-worths hide opposing groups.", ct(n_resp), resp_h)
  } else {

Uc <- matrix(0, nrow(A), length(used_attrs))
    for (j in seq_along(used_attrs)) {
      ac <- used_attrs[j]; sel <- which(lvl_attr == ac)
      Uc[, j] <- u[sel][match(A[[ac]], lvl_name[sel])]
    }
    u_total <- rowSums(Uc)
    rlist <- split(seq_along(pref), resp_vec)
    if (length(rlist) > 800) { set.seed(42); rlist <- rlist[sample(length(rlist), 800)] }
    ha <- character(0); hagg <- character(0); hn <- integer(0)
    hmatch <- integer(0); hpct <- numeric(0); hsecond <- character(0)
    for (j in seq_along(used_attrs)) {
      ac <- used_attrs[j]
      sel <- lvl_attr == ac
      uu <- u[sel]; nm <- lvl_name[sel]
      hi <- safe_which_max(uu)
      top_lv <- if (is.na(hi)) NA_character_ else nm[hi]
      partial <- pref - (u_total - Uc[, j])
      assessed <- 0L; matched <- 0L; alt_counts <- integer(0)
      for (ix in rlist) {
        lv <- A[[ac]][ix]
        if (length(unique(lv)) < 2) next
        mns <- tapply(partial[ix], lv, mean, na.rm = TRUE)
        mns <- mns[is.finite(mns)]
        if (length(mns) < 2) next
        best <- names(mns)[mns == max(mns)]
        if (length(best) != 1) next
        assessed <- assessed + 1L
        if (!is.na(top_lv) && identical(best, top_lv)) {
          matched <- matched + 1L
        } else {
          alt_counts[best] <- (if (is.na(alt_counts[best])) 0L else alt_counts[best]) + 1L
        }
      }
      if (assessed < 8) next
      ha <- c(ha, unname(attr_h[[ac]]))
      hagg <- c(hagg, top_lv %||% "n/a")
      hn <- c(hn, assessed); hmatch <- c(hmatch, matched)
      hpct <- c(hpct, 100 * matched / assessed)
      hsecond <- c(hsecond, if (length(alt_counts) == 0) "none" else
        names(alt_counts)[safe_which_max(as.numeric(alt_counts))])
    }
    if (length(ha) > 0) {
      het_ok <- TRUE
      o <- order(hpct)
      het_df <- data.frame(
        attribute = ha[o],
        agreement_pct = round(hpct[o], 1),
        aggregate_top_level = hagg[o],
        most_common_dissent = hsecond[o],
        respondents_assessed = hn[o],
        stringsAsFactors = FALSE
      )
      het_min <- het_df$agreement_pct[1]
      het_worst_attr <- het_df$attribute[1]
      het_flag <- het_min < 65
      het_df$interpretation <- ifelse(
        het_df$agreement_pct < 65,
        sprintf("Only %s%% of respondents rank &#x27;%s' top on this attribute; the most common alternative favourite is '%s'. The aggregate part-worths for this attribute average over groups that want different things.",
                vapply(het_df$agreement_pct, r2, ""), het_df$aggregate_top_level,
                het_df$most_common_dissent),
        sprintf("%s%% of respondents rank &#x27;%s' top on this attribute, so the aggregate part-worths describe most of the sample rather than blending opposed groups.",
                vapply(het_df$agreement_pct, r2, ""), het_df$aggregate_top_level))
    } else {
      het_reason <- sprintf("Each respondent in &#x27;%s' saw too few distinct levels of any attribute to work out their individual favourite, so the aggregate part-worths could not be checked for a split audience.", resp_h)
    }
  }
  if (is.null(het_df)) {
    het_df <- data.frame(
      attribute = "Heterogeneity check not run",
      agreement_pct = NA_real_, aggregate_top_level = "n/a",
      most_common_dissent = "n/a", respondents_assessed = NA_integer_,
      interpretation = het_reason, stringsAsFactors = FALSE)
  }

norm_note <- sprintf(
    "Part-worths are identified only up to a constant within each attribute, so they are reported here re-centred to sum to zero inside every attribute: each number is that level against the average level of its OWN attribute. Comparing a level of one attribute against a level of another is meaningful only through the attribute-importance ranges, never by reading the two part-worths side by side.")
  tested_note <- paste0(
    "Importance is conditional on the levels that were tested: ",
    paste(sprintf("%s across %s", imp_df$attribute, imp_df$tested_levels),
          collapse = "; "),
    ". Widen or narrow any of these ranges and the importance ordering can change — a price range of two dollars would make price look unimportant no matter how price-sensitive the market is.")
  stated_note <- paste0(
    "Every number here is a stated preference: it describes choices people made between hypothetical profiles in this survey, not purchases they made with their own money. Stated preference routinely overstates what people will actually pay and understates the friction of a real purchase.")

methods_df <- data.frame(
    item = c("Data shape detected", "Estimation", "Normalisation",
             "Attribute importance", "Tested ranges", "Willingness-to-pay",
             "Market simulation", "Heterogeneity", "What this is not"),
    detail = c(
      shape_reason,
      if (shape == "choice")
        sprintf("Conditional(multinomial) logit fitted by maximising the conditional-logit log-likelihood directly with a quasi-Newton optimiser and an analytic gradient; standard errors come from the observed information matrix. Final log-likelihood %s across %s choice tasks; %s = %s.",
                r2(ll_final), ct(n_tasks), fit_label, r3(fit_stat))
      else
        sprintf("Ordinary least squares of &#x27;%s' on dummy-coded attribute levels. %s = %s; residual standard deviation %s across %s rows.",
                pref_h, fit_label, r3(fit_stat), r3(resid_sd), ct(final_rows)),
      norm_note,
      sprintf("For each attribute, the range between its best and worst part-worth, divided by the sum of those ranges. Here the ranges are %s, summing to %s.",
              paste(sprintf("%s %s", imp_df$attribute, vapply(imp_df$utility_range, r3, "")), collapse = "; "),
              r3(tot_range)),
      tested_note,
      if (wtp_ok)
        sprintf("The price part-worths were regressed on the numeric price levels, giving a slope of %s utility per unit of &#x27;%s' (95%% interval %s to %s, linearity R-squared %s across the %s tested price points). Willingness-to-pay for a level is its part-worth divided by that slope; the interval uses the delta method with the full covariance between the level's part-worth and the slope.",
                r3(price_slope), unname(attr_h[[price_attr]]), r3(price_slope_lo),
                r3(price_slope_hi), r3(price_lin_r2), ct(length(price_levels_num)))
      else wtp_reason,
      if (shape == "choice")
        sprintf("Each configuration&#x27;s utility is the sum of its levels' part-worths; shares are the logit shares among the %s configurations listed and nothing else. They are model predictions inside the tested attribute space, not forecasts of real market share — no competitor outside this list, no distribution, no awareness, no availability.",
                ct(nrow(sim_df)))
      else
        sprintf("Each configuration&#x27;s utility is the sum of its levels' part-worths, so the configurations can be ranked and the utility gaps read. %s",
                share_basis),
      if (het_ok)
        sprintf("For each attribute, the aggregate contribution of every OTHER attribute was subtracted from &#x27;%s' row by row, and each respondent's own favourite level was read off the resulting partial scores and compared against the aggregate favourite. Agreement ranged from %s%% (%s) to %s%% across %s respondents assessed. Low agreement can also mean low power when a respondent evaluated few profiles, so read it alongside the number assessed.",
                pref_h, r2(min(het_df$agreement_pct, na.rm = TRUE)),
                het_df$attribute[1], r2(max(het_df$agreement_pct, na.rm = TRUE)),
                ct(max(het_df$respondents_assessed, na.rm = TRUE)))
      else het_reason,
      stated_note
    ),
    stringsAsFactors = FALSE
  )

metrics <- list(
    `Observations`        = final_rows,
    `Attributes`          = length(used_attrs),
    `Levels Estimated`    = length(u),
    `Data Shape`          = if (shape == "choice") "choice-based" else "ratings-based",
    `Most Important`      = sprintf("%s(%s%%)", top_imp$attribute, r2(top_imp$importance_pct)),
    `Fit`                 = sprintf("%s = %s", fit_label, r3(fit_stat))
  )
  if (shape == "choice") metrics[["Choice Tasks"]] <- as.integer(n_tasks)
  if (wtp_ok) metrics[["Value Per Utility Point"]] <- round(money_per_util, 3)
  if (het_ok) metrics[["Lowest Agreement"]] <- round(het_min, 1)

  best_row <- sim_df[1, ]
  wtp_clause <- if (wtp_ok) {
    top_w <- wtp_df[1, ]
    sprintf(" Against the price slope of %s utility per unit of &#x27;%s', the best-valued level is %s at about %s%s (95%% interval %s%s to %s%s) — a stated-preference estimate, not a price anyone has agreed to pay.",
            r3(price_slope), unname(attr_h[[price_attr]]), top_w$level_label,
            price_sym, r2(top_w$wtp), price_sym, r2(top_w$wtp_low), price_sym, r2(top_w$wtp_high))
  } else {
    paste0(" Willingness-to-pay is not reported: ", wtp_reason)
  }
  het_clause <- if (het_ok && het_flag) {
    sprintf(" Heterogeneity warning: only %s%% of respondents share the aggregate favourite level of &#x27;%s', so the aggregate part-worths for that attribute blend groups that want opposite things.",
            r2(het_min), het_worst_attr)
  } else if (het_ok) {
    sprintf(" Respondent-level agreement with the aggregate favourite runs from %s%% to %s%%, so the aggregate is not obviously masking opposed segments.",
            r2(min(het_df$agreement_pct, na.rm = TRUE)),
            r2(max(het_df$agreement_pct, na.rm = TRUE)))
  } else ""
  share_clause <- if (shape == "choice") {
    sprintf(" In a simulated choice between the %s configurations built from the tested levels, &#x27;%s' takes the largest share at %s%%.",
            ct(nrow(sim_df)), best_row$configuration, r2(best_row$predicted_share_pct))
  } else {
    sprintf(" Of the %s configurations built from the tested levels, &#x27;%s' carries the highest total utility at %s.",
            ct(nrow(sim_df)), best_row$configuration, r2(best_row$total_utility))
  }

  json_output <- list(
    answer = paste0(
      "Conjoint analysis of ", ct(final_rows), " ", if (shape == "choice")
        paste0("alternatives across ", ct(n_tasks), " choice tasks") else
        paste0("profile evaluations of &#x27;", pref_h, "'"),
      " across ", ct(length(used_attrs)), " attributes(",
      paste(unname(attr_h[used_attrs]), collapse = ", "), "): &#x27;",
      top_imp$attribute, "&#x27; is the most important attribute at ",
      r2(top_imp$importance_pct), "% of the total part-worth range, driven by the gap between &#x27;",
      top_imp$worst_level, "&#x27; and '", top_imp$best_level,
      "&#x27; — an importance that holds only for the levels actually tested (",
      top_imp$tested_levels, ").", wtp_clause, share_clause, het_clause,
      " ", stated_note
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "part_worths",
        "attribute_importance", "willingness_to_pay", "market_simulation",
        "heterogeneity", "methods"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed, n_incomplete = n_incomplete,
    bad_tasks = bad_tasks, lumped = lumped,
    shape = shape, shape_reason = shape_reason, converged = converged,
    pref_h = pref_h, task_h = task_h, resp_h = resp_h, attr_h = attr_h,
    used_attrs = used_attrs, dropped_df = dropped_df, levels_list = levels_list,
    pw_df = pw_df, imp_df = imp_df, top_imp = top_imp, tot_range = tot_range,
    wtp_df = wtp_df, wtp_ok = wtp_ok, wtp_reason = wtp_reason,
    money_per_util = money_per_util, price_sym = price_sym,
    price_attr = if (is.na(price_attr)) NA_character_ else unname(attr_h[[price_attr]]),
    price_slope = price_slope, price_slope_lo = price_slope_lo,
    price_slope_hi = price_slope_hi, price_lin_r2 = price_lin_r2,
    price_levels_num = price_levels_num,
    sim_df = sim_df, het_df = het_df, het_ok = het_ok, het_flag = het_flag,
    het_reason = het_reason,
    het_min = het_min, het_worst_attr = het_worst_attr,
    methods_df = methods_df, fit_stat = fit_stat, fit_label = fit_label,
    resid_sd = resid_sd, base_utility = base_utility, ll_final = ll_final,
    n_resp = n_resp, n_tasks = n_tasks,
    norm_note = norm_note, tested_note = tested_note, stated_note = stated_note,
    metrics = metrics, json_output = json_output
  )
}