set.seed(1)
options(digits = 5)\[ \newcommand{\x}{\mathbf{x}} \newcommand{\bfbeta}{\boldsymbol{\beta}} \newcommand{\betabar}{\bar{\bfbeta}} \newcommand{\bftheta}{\boldsymbol{\theta}} \newcommand{\iid}{\stackrel{iid}{\sim}} \newcommand{\E}{\operatorname{E}} \newcommand{\coloneqq}{\mathrel{:=}} \]
What this notebook does
The binary dual-response model is the incumbent. If a practitioner fields a five-point purchase-likelihood scale and then collapses it (“top-2-box counts as a buy”), the resulting binary model is correctly specified: collapsing a cumulative-link model at one of its own cut points gives back a cumulative-link model with one cut point. So dichotomization costs no consistency. What it costs is precision, and this notebook measures how much, three ways:
- Exactly, from the Fisher information. No simulation noise at all: the asymptotic variance of every estimator is computed in closed form across four parameter configurations.
- By paired Monte Carlo, fitting the ordinal model and every binary variant to the same simulated datasets, so nothing differs except the coarsening of the second response.
- At the individual level, in the hierarchical model a practitioner would actually estimate, where each respondent contributes only twelve tasks and pooled cut points can no longer rescue a coarsened second stage.
The third is where the argument turns. The aggregate penalties are real but modest, and an honest reading of them alone would say dichotomization is a minor sin. It is not, and the reason is that aggregate prediction variance is dominated by the forced-choice stage that both estimators share.
This notebook is the source of dichotomization_study.rds and hb_simstudy_dichot.rds.
In [1]:
1. Why coarsening is lossy but not biased
Under the model, \(\Pr(y \ge k) = 1 - G(c_{k-1} - \overline{\mu})\). Define a binary variable \(b = \mathbb{1}\{y \ge k\}\) and its likelihood is exactly a binary dual-response model with the single cut point \(c_{k-1}\). The parameters \(\bfbeta\) and \(c_{k-1}\) are still identified and still consistently estimated.
What is lost is the information in the other \(W - 2\) thresholds. The paper’s information decomposition makes this precise: the second response contributes \(\omega(\overline{\mu}; \mathbf{c}) \, \bar{x} \bar{x}'\) to the information matrix, where
\[ \omega = \sum_{w} \frac{\left( g_w - g_{w-1} \right)^2}{\pi_w} \]
sums over the categories. Dropping categories drops terms. The quantity \(\omega\) is computable in closed form, so we can say exactly what fraction of the second response’s information each coarsening retains, independent of any particular design.
2. Machinery
Copies of R/lib/dgp.R, R/lib/likelihood.R and R/homogeneous/fisher.R. The Fisher block is the part worth reading: it implements the paper’s information decomposition directly.
In [2]:
Design, simulator, likelihood, and exact Fisher information
rgumbel <- function(n) -log(-log(runif(n)))
# Attribute levels are drawn at random rather than from an efficiency-optimized
# catalogue of the kind commercial practice uses (one-way and two-way level
# balance, D-efficiency, utility balance). Nothing here needs one. The headline
# comparisons are paired -- competing specifications are fit to the same
# simulated data on the same design -- so design efficiency cancels from them,
# and recovery needs only the rank condition, which id_rank_check() verifies.
# Optimizing for the forced choice alone would not obviously serve the second
# response in any case: it is identified by variation in task attractiveness
# across tasks, which within-task utility balance does nothing to produce.
make_design <- function(n_tasks, J, seed) {
set.seed(seed)
n_rows <- n_tasks * J
a <- sample(1:3, n_rows, replace = TRUE); b <- sample(1:3, n_rows, replace = TRUE)
price <- runif(n_rows, 0.5, 2.5)
X <- cbind(a2 = as.numeric(a == 2), a3 = as.numeric(a == 3),
b2 = as.numeric(b == 2), b3 = as.numeric(b == 3), price = price)
list(X = X, n_tasks = n_tasks, J = J, P = ncol(X),
task_id = rep(seq_len(n_tasks), each = J), # task of each X row
alt_id = rep(seq_len(J), times = n_tasks)) # alternative within task
}
id_rank_check <- function(X) {
aug <- cbind(X, 1); list(rank = qr(aug)$rank, required = ncol(aug),
pass = qr(aug)$rank == ncol(aug))
}
par_to_cut <- function(par, P, W) {
if (W == 2) return(par[P + 1])
par[P + 1] + c(0, cumsum(exp(par[(P + 2):(P + W - 1)])))
}
cut_to_par <- function(cut) if (length(cut) == 1) cut else c(cut[1], log(diff(cut)))
row_max <- function(M) do.call(pmax, as.data.frame(M))
ord_prob_z <- function(lo, hi, link) {
if (link == "logit") plogis(hi) * plogis(-lo) * (-expm1(lo - hi))
else { ea <- exp(-lo); eb <- exp(-hi)
p <- exp(-eb) * (-expm1(-(ea - eb))); p[!is.finite(eb)] <- 0; p }
}
ord_prob <- function(mubar, cut, y, link) {
caug <- c(-Inf, cut, Inf); ord_prob_z(caug[y] - mubar, caug[y + 1L] - mubar, link)
}
negloglik <- function(par, dat) {
design <- dat$design; P <- design$P; W <- dat$W
n <- design$n_tasks; J <- design$J
beta <- par[1:P]; cut <- par_to_cut(par, P, W)
V <- matrix(design$X %*% beta, nrow = n, ncol = J, byrow = TRUE)
m <- row_max(V); mubar <- m + log(rowSums(exp(V - m)))
-(sum(V[cbind(seq_len(n), dat$jstar)] - mubar) +
sum(log(pmax(ord_prob(mubar, cut, dat$y, dat$link), 1e-312))))
}
num_grad <- function(f, x, eps = 1e-6) {
vapply(seq_along(x), function(k) {
h <- eps * max(1, abs(x[k])); xp <- x; xp[k] <- xp[k] + h
xm <- x; xm[k] <- xm[k] - h; (f(xp) - f(xm)) / (2 * h) }, numeric(1))
}
simulate_dual <- function(design, beta, cut, link = c("cloglog","logit"), seed) {
link <- match.arg(link); set.seed(seed)
n <- design$n_tasks; J <- design$J
V <- matrix(design$X %*% beta, nrow = n, ncol = J, byrow = TRUE)
u <- V + matrix(rgumbel(n * J), n, J)
jstar <- max.col(u, ties.method = "first")
ustar <- u[cbind(seq_len(n), jstar)]
latent <- if (link == "cloglog") ustar else ustar - rgumbel(n)
list(design = design, jstar = jstar, y = findInterval(latent, cut) + 1L,
W = length(cut) + 1L, link = link, beta_true = beta, cut_true = cut)
}
fit_dual_mle <- function(dat, start = NULL) {
design <- dat$design; P <- design$P; W <- dat$W
if (is.null(start)) {
freq <- tabulate(dat$y, nbins = W)
cumq <- pmin(pmax(cumsum(freq)[1:(W - 1)] / sum(freq), 1e-4), 1 - 1e-4)
ginv <- if (dat$link == "logit") qlogis else function(q) -log(-log(q))
cut0 <- log(design$J) + ginv(cumq)
if (W > 2) for (k in 2:(W - 1)) cut0[k] <- max(cut0[k], cut0[k - 1] + 1e-3)
start <- c(rep(0, P), cut_to_par(cut0))
}
fn <- function(p) negloglik(p, dat)
opt <- optim(start, fn, method = "BFGS", control = list(maxit = 1000, reltol = 1e-12))
H <- optimHess(opt$par, fn); ev <- eigen(H, symmetric = TRUE, only.values = TRUE)$values
cut_hat <- par_to_cut(opt$par, P, W)
Jc <- matrix(0, W - 1, length(opt$par)); Jc[, P + 1] <- 1
if (W > 2) { d <- exp(opt$par[(P + 2):(P + W - 1)])
for (w in 2:(W - 1)) Jc[w, (P + 2):(P + w)] <- d[1:(w - 1)] }
Vpar <- tryCatch(solve(H), error = function(e) matrix(NA, length(opt$par), length(opt$par)))
list(beta = opt$par[1:P], se_beta = sqrt(diag(Vpar)[1:P]), cut = cut_hat,
se_cut = sqrt(diag(Jc %*% Vpar %*% t(Jc))), nll = opt$value,
convergence = opt$convergence, grad_max = max(abs(num_grad(fn, opt$par))),
eigen = ev, min_eig = min(ev),
cond = max(ev) / max(min(ev), .Machine$double.eps), par = opt$par)
}
# ---- exact Fisher information (R/homogeneous/fisher.R) ---------------------
link_cdf <- function(z, link) if (link == "logit") plogis(z) else exp(-exp(-z))
link_pdf <- function(z, link) if (link == "logit") dlogis(z) else exp(-(z + exp(-z)))
# Information from one task: the MNL block plus the ordinal-stage block.
fisher_task <- function(Xt, beta, cut, link) {
P <- ncol(Xt); W <- length(cut) + 1L; K <- P + W - 1L
V <- as.numeric(Xt %*% beta)
eV <- exp(V - max(V)); p <- eV / sum(eV)
mubar <- max(V) + log(sum(eV))
I <- matrix(0, K, K)
I[1:P, 1:P] <- t(Xt) %*% (diag(p) - tcrossprod(p)) %*% Xt # choice stage
z <- cut - mubar
Gz <- c(0, link_cdf(z, link), 1); gz <- c(0, link_pdf(z, link), 0)
pi_w <- diff(Gz)
xbar <- as.numeric(t(p) %*% Xt)
for (w in seq_len(W)) { # ordinal stage
s <- numeric(K)
s[1:P] <- -(gz[w + 1] - gz[w]) / pi_w[w] * xbar
if (w <= W - 1) s[P + w] <- s[P + w] + gz[w + 1] / pi_w[w]
if (w >= 2) s[P + w - 1] <- s[P + w - 1] - gz[w] / pi_w[w]
I <- I + pi_w[w] * tcrossprod(s)
}
I
}
fisher_total <- function(design, beta, cut, link) {
n <- design$n_tasks; J <- design$J; K <- design$P + length(cut)
I <- matrix(0, K, K)
for (t in seq_len(n)) {
rows <- ((t - 1) * J + 1):(t * J)
I <- I + fisher_task(design$X[rows, , drop = FALSE], beta, cut, link)
}
I
}
# The scalar weight omega: how much the second response contributes per task.
omega_weight <- function(mubar, cut, link) {
z <- outer(mubar, cut, function(m, c) c - m)
Gz <- cbind(0, link_cdf(z, link), 1); gz <- cbind(0, link_pdf(z, link), 0)
W <- length(cut) + 1L
om <- 0
for (w in seq_len(W)) om <- om + (gz[, w + 1] - gz[, w])^2 / (Gz[, w + 1] - Gz[, w])
om
}In [3]:
MODEL <- "logit"; N_TASKS <- 3600; J <- 4; W <- 5; P <- 5
R_REPS <- 100; N_EVAL <- 200
MC_SETS <- c("set1_moderate", "set4_skewed")
param_sets <- list(
set1_moderate = list(beta = c(a2=0.8,a3=-0.5,b2=0.4,b3=1.0,price=-0.9),
cut = c(-1.0, 0.2, 1.2, 2.2)),
set2_negative = list(beta = c(a2=-1.2,a3=0.3,b2=-0.6,b3=0.7,price=-1.5),
cut = c(-2.5, -1.2, 0.5, 1.0)),
set3_stress = list(beta = c(a2=2.0,a3=-1.8,b2=1.2,b3=-0.7,price=-2.2),
cut = c(-3.5, -2.2, 0.8, 3.0)),
set4_skewed = list(beta = c(a2=0.8,a3=-0.5,b2=0.4,b3=1.0,price=-0.9),
cut = c(0.2, 1.2, 2.4, 3.6))
)
design <- make_design(N_TASKS, J, seed = 3001)
eval_design <- make_design(N_EVAL, J, seed = 777)
stopifnot(id_rank_check(design$X)$pass)set4_skewed is the configuration that matters most for practice. Its cut points push the response distribution downward, producing the bottom-heavy purchase intent that fielded scales actually show, and leaving the top-2-box convention badly placed.
In [4]:
mubar_of <- function(des, beta) {
V <- matrix(des$X %*% beta, nrow = des$n_tasks, ncol = des$J, byrow = TRUE)
m <- row_max(V); m + log(rowSums(exp(V - m)))
}
xbar_of <- function(des, beta) {
V <- matrix(des$X %*% beta, nrow = des$n_tasks, ncol = des$J, byrow = TRUE)
m <- row_max(V); pr <- exp(V - m) / rowSums(exp(V - m))
out <- matrix(0, des$n_tasks, ncol(des$X))
for (j in seq_len(des$J))
out <- out + pr[, j] * des$X[seq(j, by = des$J, length.out = des$n_tasks), , drop = FALSE]
out
}
exceed_pred <- function(des, beta, cut_scalar)
1 - link_cdf(cut_scalar - mubar_of(des, beta), MODEL)
# Which position in each estimator's cut vector can predict P(y >= k); NA = it cannot.
cutpos <- function(est, k) switch(est,
ord5 = k - 1, ord3 = if (k %in% 3:4) k - 2 else NA,
bin2 = if (k == 2) 1 else NA, bin3 = if (k == 3) 1 else NA,
bin4 = if (k == 4) 1 else NA, bin5 = if (k == 5) 1 else NA)
est_cut <- function(est, cut_true) switch(est,
ord5 = cut_true, ord3 = cut_true[2:3], bin2 = cut_true[1],
bin3 = cut_true[2], bin4 = cut_true[3], bin5 = cut_true[4])
EST <- c("ord5", "ord3", "bin2", "bin3", "bin4", "bin5")
grad_exceed <- function(Xt, beta, cut, pos, npar) {
V <- as.numeric(Xt %*% beta); m <- max(V)
mubar <- m + log(sum(exp(V - m))); p <- exp(V - mubar)
g <- link_pdf(cut[pos] - mubar, MODEL)
gr <- numeric(npar); gr[1:P] <- g * as.numeric(t(p) %*% Xt); gr[P + pos] <- -g
gr
}Note what cutpos encodes: a binary estimator cut at \(k\) can answer questions about threshold \(k\) and no other. A top-box model cannot tell you about top-2-box demand. This is not an efficiency loss, it is an inability, and it returns in the pricing counterfactual of notebook 08.
3. Exact efficiency, from the Fisher information
In [5]:
analytic <- list()
for (sn in names(param_sets)) {
ps <- param_sets[[sn]]
mu_est_design <- mubar_of(design, ps$beta)
catprob <- diff(c(0, colMeans(link_cdf(outer(mu_est_design, ps$cut,
function(m, c) c - m), MODEL)), 1))
u <- colMeans(xbar_of(design, ps$beta)); u <- u / sqrt(sum(u^2))
Vcov <- list(); sdtab <- list()
for (e in EST) {
ct <- est_cut(e, ps$cut)
Vc <- solve(fisher_total(design, ps$beta, ct, MODEL))
Vcov[[e]] <- Vc
Vb <- Vc[1:P, 1:P]
row <- data.frame(estimator = e, mean_sd_beta = mean(sqrt(diag(Vb))),
sd_u_beta = sqrt(t(u) %*% Vb %*% u),
mean_omega = mean(omega_weight(mu_est_design, ct, MODEL)))
for (k in 2:W) {
pos <- cutpos(e, k)
if (!is.na(pos)) {
grads <- vapply(seq_len(N_EVAL), function(t) {
rows <- ((t - 1) * J + 1):(t * J)
grad_exceed(eval_design$X[rows, , drop = FALSE], ps$beta, ct, pos, nrow(Vc))
}, numeric(nrow(Vc)))
pv <- colSums(grads * (Vc %*% grads))
gbar <- rowMeans(grads)
row[[paste0("sd_pred_k", k)]] <- sqrt(mean(pv))
row[[paste0("sd_demand_k", k)]] <- sqrt(t(gbar) %*% Vc %*% gbar)
row[[paste0("sd_c", k - 1)]] <- sqrt(Vc[P + pos, P + pos])
} else {
row[[paste0("sd_pred_k", k)]] <- NA
row[[paste0("sd_demand_k", k)]] <- NA
row[[paste0("sd_c", k - 1)]] <- NA
}
}
sdtab[[e]] <- row
}
sdtab <- do.call(rbind, sdtab); rownames(sdtab) <- NULL
analytic[[sn]] <- list(sd = sdtab, u = u, catprob = catprob)
}Implied response distributions
In [6]:
knitr::kable(
do.call(rbind, lapply(names(param_sets), function(sn) {
cp <- analytic[[sn]]$catprob
data.frame(configuration = sn, `1` = cp[1], `2` = cp[2], `3` = cp[3],
`4` = cp[4], `5` = cp[5], `top-2-box` = sum(cp[4:5]),
check.names = FALSE)})),
row.names = FALSE, digits = 3)| configuration | 1 | 2 | 3 | 4 | 5 | top-2-box |
|---|---|---|---|---|---|---|
| set1_moderate | 0.145 | 0.205 | 0.233 | 0.202 | 0.215 | 0.417 |
| set2_negative | 0.152 | 0.221 | 0.365 | 0.080 | 0.181 | 0.262 |
| set3_stress | 0.082 | 0.125 | 0.511 | 0.225 | 0.057 | 0.282 |
| set4_skewed | 0.350 | 0.233 | 0.233 | 0.119 | 0.065 | 0.184 |
What each coarsening costs
In [7]:
eff_tab <- do.call(rbind, lapply(names(param_sets), function(sn) {
a <- analytic[[sn]]$sd
data.frame(configuration = sn,
`top-2-box, beta` = round((a$mean_sd_beta[5] / a$mean_sd_beta[1])^2, 3),
`worst cut, beta` = round(max((a$mean_sd_beta / a$mean_sd_beta[1])^2, na.rm = TRUE), 3),
`omega share, top-2-box` = round(a$mean_omega[5] / a$mean_omega[1], 2),
`omega share, worst` = round(min(a$mean_omega[-1] / a$mean_omega[1], na.rm = TRUE), 2),
check.names = FALSE)}))
knitr::kable(eff_tab, row.names = FALSE)| configuration | top-2-box, beta | worst cut, beta | omega share, top-2-box | omega share, worst |
|---|---|---|---|---|
| set1_moderate | 1.041 | 1.091 | 0.73 | 0.38 |
| set2_negative | 1.088 | 1.119 | 0.59 | 0.40 |
| set3_stress | 1.152 | 1.327 | 0.59 | 0.19 |
| set4_skewed | 1.080 | 1.125 | 0.48 | 0.20 |
Two columns, two different things. The beta columns are variance ratios: a value of 1.05 means the binary analyst needs five percent more sample for the same precision on preference parameters. The omega share columns say what fraction of the second response’s own information survives the coarsening, and those numbers are much worse. The gap between them exists because \(\bfbeta\) is mostly identified by the forced choice, which both estimators keep intact.
The analyst must commit to the cut before seeing the data, so the “worst cut” column is the relevant risk, not the best case.
Full detail, per configuration
In [8]:
for (sn in names(param_sets)) {
cat("\n###", sn, "\n")
a <- analytic[[sn]]$sd
ratio <- function(col) round((a[[col]] / a[[col]][1])^2, 3)
print(knitr::kable(
data.frame(estimator = a$estimator,
beta_avg = round((a$mean_sd_beta / a$mean_sd_beta[1])^2, 3),
beta_u_dir = ratio("sd_u_beta"),
omega_share = round(a$mean_omega / a$mean_omega[1], 3),
pred_k4 = ratio("sd_pred_k4"),
demand_k4 = ratio("sd_demand_k4")),
row.names = FALSE))
}
### set1_moderate
|estimator | beta_avg| beta_u_dir| omega_share| pred_k4| demand_k4|
|:---------|--------:|----------:|-----------:|-------:|---------:|
|ord5 | 1.000| 1.000| 1.000| 1.000| 1.000|
|ord3 | 1.015| 1.014| 0.890| 1.010| 1.006|
|bin2 | 1.091| 1.085| 0.381| NA| NA|
|bin3 | 1.043| 1.040| 0.686| NA| NA|
|bin4 | 1.041| 1.038| 0.732| 1.023| 1.013|
|bin5 | 1.078| 1.074| 0.516| NA| NA|
### set2_negative
|estimator | beta_avg| beta_u_dir| omega_share| pred_k4| demand_k4|
|:---------|--------:|----------:|-----------:|-------:|---------:|
|ord5 | 1.000| 1.000| 1.000| 1.000| 1.000|
|ord3 | 1.016| 1.019| 0.920| 1.009| 1.005|
|bin2 | 1.119| 1.122| 0.397| NA| NA|
|bin3 | 1.056| 1.057| 0.698| NA| NA|
|bin4 | 1.088| 1.100| 0.589| 1.039| 1.018|
|bin5 | 1.117| 1.132| 0.460| NA| NA|
### set3_stress
|estimator | beta_avg| beta_u_dir| omega_share| pred_k4| demand_k4|
|:---------|--------:|----------:|-----------:|-------:|---------:|
|ord5 | 1.000| 1.000| 1.000| 1.000| 1.000|
|ord3 | 1.023| 1.016| 0.919| 1.020| 1.014|
|bin2 | 1.265| 1.238| 0.243| NA| NA|
|bin3 | 1.159| 1.133| 0.484| NA| NA|
|bin4 | 1.152| 1.139| 0.589| 1.073| 1.026|
|bin5 | 1.327| 1.305| 0.186| NA| NA|
### set4_skewed
|estimator | beta_avg| beta_u_dir| omega_share| pred_k4| demand_k4|
|:---------|--------:|----------:|-----------:|-------:|---------:|
|ord5 | 1.000| 1.000| 1.000| 1.000| 1.000|
|ord3 | 1.026| 1.025| 0.828| 1.008| 1.003|
|bin2 | 1.037| 1.034| 0.711| NA| NA|
|bin3 | 1.035| 1.032| 0.758| NA| NA|
|bin4 | 1.080| 1.076| 0.477| 1.023| 1.008|
|bin5 | 1.125| 1.117| 0.197| NA| NA|beta_u_dir is the variance in the “average attractiveness” direction, the rank-one direction the ordinal stage informs and the forced choice identifies worst. It is consistently the most damaged, which is exactly what the theory predicts: that is where the second response carries information the first does not.
4. Monte Carlo validation
The analytic numbers are asymptotic. Do they describe finite samples? Fit every estimator to the same simulated datasets and compare empirical spread to the predicted spread.
In [9]:
collapse3 <- function(y) 1L + (y >= 3L) + (y >= 4L)
mc <- list()
for (sn in MC_SETS) {
ps <- param_sets[[sn]]
mubar_eval_true <- mubar_of(eval_design, ps$beta)
true_ex <- sapply(2:W, function(k) 1 - link_cdf(ps$cut[k - 1] - mubar_eval_true, MODEL))
beta_est <- array(NA_real_, c(R_REPS, P, length(EST)),
dimnames = list(NULL, names(ps$beta), EST))
c_own <- array(NA_real_, c(R_REPS, 4, length(EST)),
dimnames = list(NULL, paste0("c", 1:4), EST))
demand <- array(NA_real_, c(R_REPS, W - 1, length(EST)),
dimnames = list(NULL, paste0("k", 2:W), EST))
pred_sse <- matrix(0, length(EST), W - 1, dimnames = list(EST, paste0("k", 2:W)))
pred_n <- pred_sse
for (r in seq_len(R_REPS)) {
dat5 <- simulate_dual(design, ps$beta, ps$cut, link = MODEL, seed = 60000 + r)
fits <- list()
fits$ord5 <- fit_dual_mle(dat5)
dat3 <- dat5; dat3$y <- collapse3(dat5$y); dat3$W <- 3L
fits$ord3 <- fit_dual_mle(dat3)
for (k in 2:W) {
datb <- dat5; datb$y <- 1L + (dat5$y >= k); datb$W <- 2L
fits[[paste0("bin", k)]] <- fit_dual_mle(datb)
}
for (e in EST) {
beta_est[r, , e] <- fits[[e]]$beta
for (k in 2:W) {
pos <- cutpos(e, k); if (is.na(pos)) next
c_own[r, k - 1, e] <- fits[[e]]$cut[pos]
pr <- exceed_pred(eval_design, fits[[e]]$beta, fits[[e]]$cut[pos])
demand[r, k - 1, e] <- mean(pr)
pred_sse[e, k - 1] <- pred_sse[e, k - 1] + mean((pr - true_ex[, k - 1])^2)
pred_n[e, k - 1] <- pred_n[e, k - 1] + 1
}
}
}
pred_rmse <- sqrt(pred_sse / pmax(pred_n, 1)); pred_rmse[pred_n == 0] <- NA
vr <- sapply(2:W, function(k)
(pred_rmse[paste0("bin", k), k - 1] / pred_rmse["ord5", k - 1])^2)
names(vr) <- paste0("k", 2:W)
mc[[sn]] <- list(beta_est = beta_est, c_own = c_own, demand = demand,
pred_rmse = pred_rmse, pred_var_ratio = vr, true_ex = true_ex)
}Empirical versus analytic
In [10]:
for (sn in MC_SETS) {
cat("\n###", sn, "-- cut-point sd, empirical vs analytic\n")
rows <- list()
for (k in 2:W) for (e in EST) {
if (!is.na(cutpos(e, k))) {
rows[[length(rows) + 1]] <- data.frame(
cut = paste0("c", k - 1), estimator = e,
empirical = round(sd(mc[[sn]]$c_own[, k - 1, e]), 4),
analytic = round(analytic[[sn]]$sd[[paste0("sd_c", k - 1)]][match(e, EST)], 4))
}
}
print(knitr::kable(do.call(rbind, rows), row.names = FALSE))
}
### set1_moderate -- cut-point sd, empirical vs analytic
|cut |estimator | empirical| analytic|
|:---|:---------|---------:|--------:|
|c1 |ord5 | 0.0749| 0.0783|
|c1 |bin2 | 0.0748| 0.0806|
|c2 |ord5 | 0.0669| 0.0713|
|c2 |ord3 | 0.0694| 0.0717|
|c2 |bin3 | 0.0673| 0.0725|
|c3 |ord5 | 0.0669| 0.0705|
|c3 |ord3 | 0.0692| 0.0710|
|c3 |bin4 | 0.0720| 0.0716|
|c4 |ord5 | 0.0747| 0.0740|
|c4 |bin5 | 0.0765| 0.0759|
### set4_skewed -- cut-point sd, empirical vs analytic
|cut |estimator | empirical| analytic|
|:---|:---------|---------:|--------:|
|c1 |ord5 | 0.0686| 0.0714|
|c1 |bin2 | 0.0673| 0.0725|
|c2 |ord5 | 0.0683| 0.0707|
|c2 |ord3 | 0.0714| 0.0714|
|c2 |bin3 | 0.0720| 0.0716|
|c3 |ord5 | 0.0757| 0.0754|
|c3 |ord3 | 0.0784| 0.0761|
|c3 |bin4 | 0.0774| 0.0774|
|c4 |ord5 | 0.0964| 0.0916|
|c4 |bin5 | 0.0984| 0.0941|The two columns track each other closely, which validates both the Fisher calculation and the standard-error formulas.
Is anything biased?
In [11]:
bias_tab <- do.call(rbind, lapply(MC_SETS, function(sn) {
ps <- param_sets[[sn]]
data.frame(configuration = sn,
t(sapply(EST, function(e) {
b <- colMeans(mc[[sn]]$beta_est[, , e]) - ps$beta
max(abs(b / (apply(mc[[sn]]$beta_est[, , e], 2, sd) / sqrt(R_REPS))))
})))}))
knitr::kable(bias_tab, row.names = FALSE, digits = 2)| configuration | ord5 | ord3 | bin2 | bin3 | bin4 | bin5 |
|---|---|---|---|---|---|---|
| set1_moderate | 1.12 | 1.21 | 1.17 | 1.35 | 1.22 | 1.34 |
| set4_skewed | 1.05 | 1.14 | 1.35 | 1.22 | 1.17 | 1.04 |
Maximum \(|t|\) statistics for bias, by estimator. Nothing systematic anywhere, confirming that dichotomization costs efficiency and not consistency. This is the claim that lets the rest of the comparison be about precision alone.
5. The individual level, where it matters
Everything above pools 3,600 tasks and estimates one parameter vector. That is not what practitioners do. They estimate a hierarchical model in which each respondent contributes twelve tasks, and their deliverable is individual-level part-worths and individual purchase probabilities.
At twelve tasks per person, there is no pooling to rescue the coarsened second stage.
In [12]:
Panel simulation and the hierarchical sampler
make_panel_design <- function(N, T_tasks, J, seed) {
des <- make_design(N * T_tasks, J, seed = seed)
des$N <- N; des$T_tasks <- T_tasks
des$resp <- rep(seq_len(N), each = T_tasks)
des$respondent_id <- rep(des$resp, each = J)
des
}
draw_betas <- function(N, beta_bar, Sigma, seed) {
set.seed(seed); P <- length(beta_bar)
sweep(matrix(rnorm(N * P), N, P) %*% chol(Sigma), 2, beta_bar, "+")
}
simulate_dual_hb <- function(design, Bmat, cut, link = c("cloglog","logit"), seed) {
link <- match.arg(link); set.seed(seed)
n <- design$n_tasks; J <- design$J; N <- design$N
cutmat <- if (is.matrix(cut)) cut else matrix(cut, N, length(cut), byrow = TRUE)
W <- ncol(cutmat) + 1L
Bx <- Bmat[design$respondent_id, , drop = FALSE]
V <- matrix(rowSums(design$X * Bx), nrow = n, ncol = J, byrow = TRUE)
u <- V + matrix(rgumbel(n * J), n, J)
jstar <- max.col(u, ties.method = "first")
ustar <- u[cbind(seq_len(n), jstar)]
latent <- if (link == "cloglog") ustar else ustar - rgumbel(n)
y <- 1L + rowSums(latent >= cutmat[design$resp, , drop = FALSE])
list(design = design, resp = design$resp, respondent_id = design$respondent_id, N = N,
jstar = jstar, y = as.integer(y), W = W, link = link,
Bmat_true = Bmat, cut_true = cut)
}
split_holdout <- function(dat, n_holdout) {
des <- dat$design; T_tasks <- des$T_tasks
keep_t <- rep(seq_len(T_tasks) <= T_tasks - n_holdout, times = des$N)
sub <- function(keep) {
rows <- rep(keep, each = des$J); d <- des
d$X <- des$X[rows, , drop = FALSE]; d$n_tasks <- sum(keep)
d$T_tasks <- d$n_tasks / des$N; d$resp <- des$resp[keep]
d$respondent_id <- rep(d$resp, each = des$J)
out <- dat; out$design <- d; out$resp <- d$resp; out$respondent_id <- d$respondent_id
out$jstar <- dat$jstar[keep]; out$y <- dat$y[keep]; out
}
list(est = sub(keep_t), holdout = sub(!keep_t))
}
phi_to_cutmat <- function(Phi, P, W) {
c1 <- Phi[, P + 1]
if (W == 2) return(matrix(c1, ncol = 1))
D <- exp(Phi[, (P + 2):(P + W - 1), drop = FALSE])
cm <- matrix(0, nrow(Phi), W - 1); cm[, 1] <- c1
for (w in 2:(W - 1)) cm[, w] <- cm[, w - 1] + D[, w - 1]
cm
}
resp_loglik_all <- function(dat, Bmat, cutmat, choice_only = FALSE) {
des <- dat$design; n <- des$n_tasks; J <- des$J
Bx <- Bmat[dat$respondent_id, , drop = FALSE]
V <- matrix(rowSums(des$X * Bx), nrow = n, ncol = J, byrow = TRUE)
m <- row_max(V); mubar <- m + log(rowSums(exp(V - m)))
ll <- V[cbind(seq_len(n), dat$jstar)] - mubar
if (!choice_only) {
ca <- cbind(-Inf, cutmat, Inf)
ll <- ll + log(pmax(ord_prob_z(ca[cbind(dat$resp, dat$y)] - mubar,
ca[cbind(dat$resp, dat$y + 1L)] - mubar,
dat$link), 1e-312))
}
as.numeric(rowsum(ll, dat$resp, reorder = TRUE))
}
default_hb_priors <- function(d) list(phibar0 = rep(0, d), A = 1/100,
nu = d + 3, V0 = (d + 3) * diag(d), zeta_prec = 1/100)
draw_phibar <- function(Phi, Sigma, pr) {
N <- nrow(Phi); Si <- chol2inv(chol(Sigma))
Vb <- chol2inv(chol(N * Si + pr$A * diag(ncol(Phi))))
m <- Vb %*% (Si %*% (N * colMeans(Phi)) + pr$A * pr$phibar0)
as.numeric(m + t(chol(Vb)) %*% rnorm(ncol(Phi)))
}
draw_sigma <- function(Phi, phibar, pr) {
Vn <- pr$V0 + crossprod(sweep(Phi, 2, phibar))
chol2inv(chol(rWishart(1, pr$nu + nrow(Phi), chol2inv(chol(Vn)))[,,1]))
}
fit_dual_hb <- function(dat, mcmc = list(R=30000, burn=10000, thin=10),
het_cut = FALSE, choice_only = FALSE, priors = NULL,
keep_phi = TRUE, seed = NULL, verbose = TRUE) {
if (!is.null(seed)) set.seed(seed)
des <- dat$design; N <- dat$N; P <- des$P; W <- dat$W
dim_phi <- if (het_cut) P + W - 1 else P
if (is.null(priors)) priors <- default_hb_priors(dim_phi)
R_iter <- mcmc$R; burn <- mcmc$burn; thin <- mcmc$thin
nkeep <- floor((R_iter - burn) / thin)
Phi <- matrix(0, N, dim_phi)
freq <- tabulate(dat$y, nbins = W)
cumq <- pmin(pmax(cumsum(freq)[1:(W-1)]/sum(freq), 1e-4), 1-1e-4)
ginv <- if (dat$link == "logit") qlogis else function(q) -log(-log(q))
cut0 <- log(des$J) + ginv(cumq)
if (W > 2) for (k in 2:(W-1)) cut0[k] <- max(cut0[k], cut0[k-1] + 1e-3)
zeta <- cut_to_par(cut0)
if (het_cut) Phi[, (P+1):(P+W-1)] <- matrix(zeta, N, W-1, byrow = TRUE)
phibar <- colMeans(Phi); Sigma <- diag(dim_phi)
cur_cutmat <- if (het_cut) phi_to_cutmat(Phi, P, W) else
matrix(par_to_cut(zeta, 0, W), N, W-1, byrow = TRUE)
cur_ll <- resp_loglik_all(dat, Phi[, 1:P, drop=FALSE], cur_cutmat, choice_only)
s_phi <- rep(2.93/sqrt(dim_phi), N); acc_phi <- rep(0, N)
s_zeta <- 0.05; acc_zeta <- 0; window <- 100
keep_betabar <- matrix(NA_real_, nkeep, dim_phi)
keep_Sigma <- array(NA_real_, c(nkeep, dim_phi, dim_phi))
keep_cut <- if (!het_cut && !choice_only) matrix(NA_real_, nkeep, W-1) else NULL
keep_ll <- numeric(nkeep); keep_ll_resp <- matrix(NA_real_, nkeep, N)
keep_Phi <- if (keep_phi) array(NA_real_, c(nkeep, N, dim_phi)) else NULL
t0 <- Sys.time(); ki <- 0
for (it in seq_len(R_iter)) {
U <- chol(Sigma)
Prop <- Phi + (matrix(rnorm(N*dim_phi), N, dim_phi) %*% U) * s_phi
prop_cutmat <- if (het_cut) phi_to_cutmat(Prop, P, W) else cur_cutmat
prop_ll <- resp_loglik_all(dat, Prop[, 1:P, drop=FALSE], prop_cutmat, choice_only)
Si <- chol2inv(chol(Sigma))
dcur <- sweep(Phi, 2, phibar); dprop <- sweep(Prop, 2, phibar)
la <- (prop_ll - 0.5*rowSums((dprop %*% Si)*dprop)) -
(cur_ll - 0.5*rowSums((dcur %*% Si)*dcur))
acc <- log(runif(N)) < la
Phi[acc, ] <- Prop[acc, ]; cur_ll[acc] <- prop_ll[acc]
if (het_cut && any(acc)) cur_cutmat[acc, ] <- prop_cutmat[acc, , drop=FALSE]
acc_phi <- acc_phi + acc
if (!het_cut && !choice_only) {
zp <- zeta + rnorm(W-1, sd = s_zeta)
cp <- matrix(par_to_cut(zp, 0, W), N, W-1, byrow = TRUE)
lp_ll <- resp_loglik_all(dat, Phi[, 1:P, drop=FALSE], cp, choice_only)
if (log(runif(1)) < (sum(lp_ll) - 0.5*priors$zeta_prec*sum(zp^2)) -
(sum(cur_ll) - 0.5*priors$zeta_prec*sum(zeta^2))) {
zeta <- zp; cur_cutmat <- cp; cur_ll <- lp_ll; acc_zeta <- acc_zeta + 1
}
}
phibar <- draw_phibar(Phi, Sigma, priors)
Sigma <- draw_sigma(Phi, phibar, priors)
if (it <= burn && it %% window == 0) {
s_phi <- pmin(pmax(s_phi*exp(1.5*(acc_phi/window - 0.23)), 0.01), 10)
acc_phi <- rep(0, N)
if (!het_cut && !choice_only) {
s_zeta <- min(max(s_zeta*exp(1.5*(acc_zeta/window - 0.30)), 1e-3), 2)
acc_zeta <- 0
}
}
if (it == burn) { acc_phi <- rep(0, N); acc_zeta <- 0 }
if (it > burn && (it - burn) %% thin == 0) {
ki <- ki + 1
keep_betabar[ki, ] <- phibar; keep_Sigma[ki, , ] <- Sigma
if (!is.null(keep_cut)) keep_cut[ki, ] <- par_to_cut(zeta, 0, W)
keep_ll[ki] <- sum(cur_ll); keep_ll_resp[ki, ] <- cur_ll
if (keep_phi) keep_Phi[ki, , ] <- Phi
}
}
post <- R_iter - burn
structure(list(draws = list(phibar = keep_betabar, Sigma = keep_Sigma,
cut = keep_cut, Phi = keep_Phi),
ll = keep_ll, ll_resp = keep_ll_resp,
accept = list(phi = mean(acc_phi/post),
zeta = if (het_cut || choice_only) NA else acc_zeta/post),
settings = list(mcmc = mcmc, het_cut = het_cut,
choice_only = choice_only, P = P, W = W, N = N,
link = dat$link, priors = priors),
runtime_min = as.numeric(difftime(Sys.time(), t0, units="mins"))),
class = "dr_hb_fit")
}
hb_beta_i <- function(fit)
apply(fit$draws$Phi[, , 1:fit$settings$P, drop = FALSE], c(2,3), mean)In [13]:
N <- 300; T_EST <- 12; T_HOLD <- 2
beta_bar <- c(a2=0.8, a3=-0.5, b2=0.4, b3=1.0, price=-0.9)
Sigma_true <- diag(c(0.6, 0.6, 0.6, 0.6, 0.3))
cut_true <- c(-1.0, 0.2, 1.2, 2.2)
KBUY <- 4 # "purchase" means top-2-box throughout
REPS <- 20
ITER <- 20000
MCMC <- list(R = ITER, burn = floor(ITER * 0.4), thin = 10)
gen_panel <- function(seed, link = "logit") {
des <- make_panel_design(N, T_EST + T_HOLD, J, seed = seed)
B <- draw_betas(N, beta_bar, Sigma_true, seed = seed + 1)
split_holdout(simulate_dual_hb(des, B, cut_true, link = link, seed = seed + 2), T_HOLD)
}
task_mubar <- function(design, Bmat) {
Bx <- Bmat[design$respondent_id, , drop = FALSE]
V <- matrix(rowSums(design$X * Bx), ncol = design$J, byrow = TRUE)
m <- row_max(V); m + log(rowSums(exp(V - m)))
}
true_exceed_panel <- function(sp) {
d <- sp$holdout$design
mubar <- task_mubar(d, sp$holdout$Bmat_true)
cutm <- matrix(sp$holdout$cut_true, N, W - 1, byrow = TRUE)
1 - plogis(cutm[d$resp, KBUY - 1] - mubar)
}
pred_exceed_panel <- function(fit, sp, k = KBUY) {
d <- sp$holdout$design
mubar <- task_mubar(d, hb_beta_i(fit))
cutvec <- colMeans(fit$draws$cut)
cc <- rep(cutvec[min(k - 1, length(cutvec))], length(mubar))
G <- if (fit$settings$link == "logit") plogis else function(z) exp(-exp(-z))
1 - G(cc - mubar)
}
hit_rate <- function(fit, sp) {
d <- sp$holdout$design
Bx <- hb_beta_i(fit)[d$respondent_id, , drop = FALSE]
V <- matrix(rowSums(d$X * Bx), ncol = d$J, byrow = TRUE)
mean(max.col(V, ties.method = "first") == sp$holdout$jstar)
}
dichotomize <- function(dat, k = KBUY) { dat$y <- 1L + (dat$y >= k); dat$W <- 2L; dat }Twenty paired replications. Each generates one panel, then fits the ordinal model and the top-2-box binary model to the same forced-choice data, so every difference is attributable to the coarsened second response.
In [14]:
m <- matrix(NA_real_, REPS, 8, dimnames = list(NULL, c(
"rmse_bi_ord","rmse_bi_bin","hit_ord","hit_bin",
"rmse_ex_ord","rmse_ex_bin","sd_c3_ord","sd_c3_bin")))
for (r in seq_len(REPS)) {
sp <- gen_panel(seed = 60000 + 10 * r)
f_ord <- fit_dual_hb(sp$est, mcmc = MCMC, seed = 61000 + r, verbose = FALSE)
f_bin <- fit_dual_hb(dichotomize(sp$est), mcmc = MCMC, seed = 62000 + r, verbose = FALSE)
tex <- true_exceed_panel(sp)
m[r, ] <- c(
sqrt(mean((hb_beta_i(f_ord) - sp$est$Bmat_true)^2)),
sqrt(mean((hb_beta_i(f_bin) - sp$est$Bmat_true)^2)),
hit_rate(f_ord, sp), hit_rate(f_bin, sp),
sqrt(mean((pred_exceed_panel(f_ord, sp) - tex)^2)),
sqrt(mean((pred_exceed_panel(f_bin, sp, k = 2) - tex)^2)),
sd(f_ord$draws$cut[, KBUY - 1]), sd(f_bin$draws$cut[, 1]))
}
res_dichot <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS))
knitr::kable(rbind(mean = res_dichot$mean, mc_se = res_dichot$se), digits = 4)| rmse_bi_ord | rmse_bi_bin | hit_ord | hit_bin | rmse_ex_ord | rmse_ex_bin | sd_c3_ord | sd_c3_bin | |
|---|---|---|---|---|---|---|---|---|
| mean | 0.5237 | 0.5345 | 0.5106 | 0.5099 | 0.1114 | 0.1214 | 0.0763 | 0.0797 |
| mc_se | 0.0036 | 0.0035 | 0.0048 | 0.0045 | 0.0011 | 0.0012 | 0.0009 | 0.0007 |
In [15]:
rm_ <- res_dichot$mean
knitr::kable(data.frame(
quantity = c("individual part-worth RMSE", "individual purchase-probability RMSE",
"posterior sd of the shared cut point", "holdout hit rate"),
ordinal = c(rm_["rmse_bi_ord"], rm_["rmse_ex_ord"], rm_["sd_c3_ord"], rm_["hit_ord"]),
binary = c(rm_["rmse_bi_bin"], rm_["rmse_ex_bin"], rm_["sd_c3_bin"], rm_["hit_bin"]),
`variance ratio` = c((rm_["rmse_bi_bin"]/rm_["rmse_bi_ord"])^2,
(rm_["rmse_ex_bin"]/rm_["rmse_ex_ord"])^2,
(rm_["sd_c3_bin"]/rm_["sd_c3_ord"])^2, NA),
check.names = FALSE), row.names = FALSE, digits = 4)| quantity | ordinal | binary | variance ratio |
|---|---|---|---|
| individual part-worth RMSE | 0.5237 | 0.5345 | 1.0417 |
| individual purchase-probability RMSE | 0.1114 | 0.1214 | 1.1884 |
| posterior sd of the shared cut point | 0.0763 | 0.0797 | 1.0904 |
| holdout hit rate | 0.5106 | 0.5099 | NA |
Read the last column. The cut point and the individual purchase probabilities, which are what the dual response exists to measure, degrade substantially. The holdout hit rate does not move at all, and should not: the first response carries the ranking information and dichotomization leaves it untouched.
That is the finding. Dichotomization does no harm to what the forced choice already tells you, and real harm to the one thing the second question was added to measure.
6. Results written
In [16]:
# Quarto executes some passes from the project root and others from this
# file's directory, so anchor the output path on the project marker rather
# than trusting the working directory.
PROJ <- if (file.exists("_quarto.yml")) "." else ".."
OUT <- file.path(PROJ, "R", "output")
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)
saveRDS(list(analytic = analytic, mc = mc,
settings = list(MODEL = MODEL, N_TASKS = N_TASKS, J = J, W = W,
R_REPS = R_REPS, param_sets = param_sets)),
file.path(OUT, "dichotomization_study.rds"))
saveRDS(list(cell = "dichot", res = res_dichot,
settings = list(N = N, T_EST = T_EST, T_HOLD = T_HOLD, J = J, W = W,
REPS = REPS, ITER = ITER, KBUY = KBUY)),
file.path(OUT, "hb_simstudy_dichot.rds"))
cat("wrote dichotomization_study.rds and hb_simstudy_dichot.rds\n")wrote dichotomization_study.rds and hb_simstudy_dichot.rds