--- title: "Misspecification Robustness" subtitle: "Five ways the model could be wrong, and what each one costs" author: - name: Dan Yavorsky - name: Geoffery Zheng date: last-modified format: html: toc: true toc-depth: 3 code-fold: false execute: warning: false message: false --- {{< include ../latex_commands.qmd >}} ## What this notebook does Every result so far assumes the model is true. This notebook breaks that assumption five times, in the five ways a sceptical reader would break it, and measures the damage. In each case the data are generated from something the model does not describe, and then fit with the model anyway. | # | the data really are... | the worry | |---|---|---| | 1 | generated by the other link | we cannot tell the links apart, so the behavioral story is untestable | | 2 | generated with respondent-specific cut points | people use scales differently and we pool them | | 3 | generated with normal rather than Gumbel errors | the whole derivation rests on Gumbel | | 4 | generated with correlated alternatives | independence of irrelevant alternatives is false in real choice sets | | 5 | generated with the outside good drifting across tasks | respondents change as they work through a survey | A recurring pattern emerges, and it is worth stating before the evidence rather than after: **the damage lands on the cut points, not on the preference parameters or the demand predictions built from them.** That is reassuring for the quantities practitioners use and a genuine caveat for the structural reading of the scale labels. This notebook is the source of five results files, `hb_simstudy_{wronglink,hetcut,nongumbel,nested,drift}.rds`. ```{r} #| label: setup set.seed(1) options(digits = 5) ``` ## 1. Machinery The simulator here is the fullest version: it carries normal errors, nested (GEV) errors, heterogeneous cut points, and two kinds of task-order drift. The sampler gains the heterogeneous-cut-point mode. ```{r} #| label: machinery #| code-fold: true #| code-summary: "Designs, GEV errors, the full panel simulator, and the sampler" rgumbel <- function(n) -log(-log(runif(n))) # Positive-stable draws (Chambers-Mallows-Stuck). Conditional on a positive # stable draw, within-nest errors are Gumbel with a common shifted location, # which reproduces the nested-logit joint distribution while leaving each # marginal standard Gumbel -- so nested and independent errors are on the same # scale and need no variance matching. rposstable <- function(n, alpha) { U <- runif(n, 0, pi); W <- rexp(n) sin(alpha * U) / (sin(U))^(1 / alpha) * (sin((1 - alpha) * U) / W)^((1 - alpha) / alpha) } rnested <- function(n, J, nest) { eta <- matrix(NA_real_, n, J) for (g in sort(unique(nest$groups))) { idx <- which(nest$groups == g); lam <- nest$lambda[g] loc <- lam * log(rposstable(n, lam)) for (j in idx) eta[, j] <- loc + lam * rgumbel(n) } eta } # 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 } 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 } 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 } } 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, shock = c("gumbel","normal"), drift = NULL, nest = NULL) { link <- match.arg(link); shock <- match.arg(shock) 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 <- if (!is.null(nest)) V + rnested(n, J, nest) else { eta <- if (shock == "gumbel") rgumbel(n * J) else rnorm(n * J, sd = pi / sqrt(6)) V + matrix(eta, n, J) } jstar <- max.col(u, ties.method = "first") ustar <- u[cbind(seq_len(n), jstar)] m <- row_max(V); mubar <- m + log(rowSums(exp(V - m))) latent <- if (link == "cloglog") ustar else ustar - rgumbel(n) shift <- numeric(n) if (!is.null(drift)) { T_tasks <- design$T_tasks tt <- rep(seq_len(T_tasks), times = N) if (drift$type == "linear") shift <- drift$delta * (tt - 1) else if (drift$type == "update") { psi <- numeric(n) for (i in seq_len(N)) { rows <- which(design$resp == i) for (k in seq_along(rows)[-1]) psi[rows[k]] <- if (runif(1) < plogis(drift$rho)) mubar[rows[k-1]] else psi[rows[k-1]] } shift <- psi } } eff_cut <- cutmat[design$resp, , drop = FALSE] + shift list(design = design, resp = design$resp, respondent_id = design$respondent_id, N = N, jstar = jstar, y = as.integer(1L + rowSums(latent >= eff_cut)), W = W, link = link, Bmat_true = Bmat, cut_true = cut, shock = shock, drift = drift, nest = nest) } 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)))) as.numeric(Vb %*% (Si %*% (N * colMeans(Phi)) + pr$A * pr$phibar0) + 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) logmeanexp <- function(x) { m <- max(x); m + log(mean(exp(x - m))) } waic <- function(ll) { lppd <- sum(apply(ll, 2, logmeanexp)); p <- sum(apply(ll, 2, var)) list(waic = -2*(lppd - p), lppd = lppd, p_waic = p, elpd = lppd - p) } ``` ```{r} #| label: config N <- 300; T_EST <- 12; T_HOLD <- 2; J <- 4; W <- 5 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) P <- length(beta_bar) KBUY <- 4 # "purchase" means top-2-box throughout REPS <- 12 ITER <- 20000 MCMC <- list(R = ITER, burn = floor(ITER * 0.4), thin = 10) gen_panel <- function(seed, link = "logit", cutmat = NULL, ...) { des <- make_panel_design(N, T_EST + T_HOLD, J, seed = seed) B <- draw_betas(N, beta_bar, Sigma_true, seed = seed + 1) cutarg <- if (is.null(cutmat)) cut_true else cutmat split_holdout(simulate_dual_hb(des, B, cutarg, 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))) } pred_exceed_panel <- function(fit, sp, k = KBUY) { d <- sp$holdout$design mubar <- task_mubar(d, hb_beta_i(fit)) if (fit$settings$het_cut) { cutm <- phi_to_cutmat(apply(fit$draws$Phi, c(2,3), mean), fit$settings$P, fit$settings$W) cc <- cutm[d$resp, k - 1] } else { cv <- colMeans(fit$draws$cut); cc <- rep(cv[min(k - 1, length(cv))], 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) } res_all <- list() ``` ## 2. Wrong link Fit the logistic link to data from our model and vice versa. The question is whether the two links are distinguishable, and whether getting it wrong matters for demand. ```{r} #| label: wronglink m <- matrix(NA_real_, REPS, 6, dimnames = list(NULL, c( "bias_ex_correct_cloglog","bias_ex_wrong_cloglog","bias_ex_correct_logit","bias_ex_wrong_logit", "waic_picks_cloglog","waic_picks_logit"))) for (r in seq_len(REPS)) { row <- numeric(6) for (truth in c("cloglog","logit")) { sp <- gen_panel(seed = 70000 + 20*r + (truth == "cloglog"), link = truth) d <- sp$holdout$design mu_true <- task_mubar(d, sp$holdout$Bmat_true) Gt <- if (truth == "logit") plogis else function(z) exp(-exp(-z)) tex <- 1 - Gt(cut_true[KBUY - 1] - mu_true) fits <- list() for (fl in c("cloglog","logit")) { datf <- sp$est; datf$link <- fl fits[[fl]] <- fit_dual_hb(datf, mcmc = MCMC, seed = 71000 + 2*r + (fl == "cloglog"), verbose = FALSE) } wrong <- setdiff(c("cloglog","logit"), truth) bc <- mean(pred_exceed_panel(fits[[truth]], sp) - tex) bw <- mean(pred_exceed_panel(fits[[wrong]], sp) - tex) picked <- if (waic(fits$cloglog$ll_resp)$elpd > waic(fits$logit$ll_resp)$elpd) "cloglog" else "logit" if (truth == "cloglog") { row[1:2] <- c(bc, bw); row[5] <- picked == "cloglog" } else { row[3:4] <- c(bc, bw); row[6] <- picked == "logit" } } m[r, ] <- row } res_all$wronglink <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS)) knitr::kable(rbind(mean = res_all$wronglink$mean, mc_se = res_all$wronglink$se), digits = 4) ``` ```{r} #| label: wronglink-read w <- res_all$wronglink$mean cat(sprintf("demand bias, correct link: up to %.4f on the probability scale\n", max(abs(c(w["bias_ex_correct_cloglog"], w["bias_ex_correct_logit"]))))) cat(sprintf("demand bias, wrong link : up to %.4f\n", max(abs(c(w["bias_ex_wrong_cloglog"], w["bias_ex_wrong_logit"]))))) cat(sprintf("WAIC selects the true link in %.0f%% of replications\n", 100 * mean(c(w["waic_picks_cloglog"], w["waic_picks_logit"])))) ``` Link misspecification is a second-order concern *for demand*. The two links agree closely over the range of inclusive values this design produces, and diverge only in the tails, which is where extrapolation lives. Notebook 07 shows that the ability to tell the links apart is governed by how much the design varies task attractiveness, not by sample size. ## 3. Ignored scale-usage heterogeneity Now respondents genuinely differ in how they use the scale: each gets her own cut points. Fit the common-cut model (wrong) and the heterogeneous-cut model (right), and compare. ```{r} #| label: hetcut m <- matrix(NA_real_, REPS, 6, dimnames = list(NULL, c( "rmse_bi_common","rmse_bi_het","rmse_ex_common","rmse_ex_het", "hit_common","hit_het"))) theta_mu <- c(-1, log(1.2), log(1.0), log(1.0)) theta_sd <- c(0.35, 0.25, 0.25, 0.25) for (r in seq_len(REPS)) { set.seed(80000 + 10*r) Theta <- sweep(matrix(rnorm(N*4), N, 4) %*% diag(theta_sd), 2, theta_mu, "+") cutmat <- t(apply(Theta, 1, function(th) th[1] + c(0, cumsum(exp(th[2:4]))))) sp <- gen_panel(seed = 80000 + 10*r + 1, cutmat = cutmat) f_com <- fit_dual_hb(sp$est, mcmc = MCMC, seed = 81000 + r, verbose = FALSE) f_het <- fit_dual_hb(sp$est, mcmc = MCMC, het_cut = TRUE, seed = 82000 + r, verbose = FALSE) d <- sp$holdout$design tex <- 1 - plogis(cutmat[d$resp, KBUY - 1] - task_mubar(d, sp$holdout$Bmat_true)) m[r, ] <- c(sqrt(mean((hb_beta_i(f_com) - sp$est$Bmat_true)^2)), sqrt(mean((hb_beta_i(f_het) - sp$est$Bmat_true)^2)), sqrt(mean((pred_exceed_panel(f_com, sp) - tex)^2)), sqrt(mean((pred_exceed_panel(f_het, sp) - tex)^2)), hit_rate(f_com, sp), hit_rate(f_het, sp)) } res_all$hetcut <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS), theta_mu = theta_mu, theta_sd = theta_sd) knitr::kable(rbind(mean = res_all$hetcut$mean, mc_se = res_all$hetcut$se), digits = 4) ``` The common-cut model loses almost nothing, even though the data-generating process really does give every respondent her own thresholds with a standard deviation of 0.35 on the first one. At twelve tasks per person there is not enough information to estimate four extra parameters per respondent, so the extra flexibility buys little and costs variance. The heterogeneous-cut extension is robustness insurance and a diagnostic, not a default. ## 4. Non-Gumbel errors The entire derivation rests on Gumbel errors. Replace them with normals matched on variance and refit. ```{r} #| label: nongumbel m <- matrix(NA_real_, REPS, 4, dimnames = list(NULL, c("rmse_bi","rmse_ex","hit","bias_c3"))) for (r in seq_len(REPS)) { sp <- gen_panel(seed = 90000 + 10*r, shock = "normal") f <- fit_dual_hb(sp$est, mcmc = MCMC, seed = 91000 + r, verbose = FALSE) # With normal errors there is no closed form for demand, so the truth is # obtained by simulating the behavioral process at the holdout tasks. set.seed(92000 + r) d2 <- sp$holdout$design Bx <- sp$holdout$Bmat_true[d2$respondent_id, , drop = FALSE] V <- matrix(rowSums(d2$X * Bx), ncol = J, byrow = TRUE) cnt <- numeric(nrow(V)) for (mc in seq_len(1000)) { u <- V + matrix(rnorm(length(V), sd = pi/sqrt(6)), nrow(V), J) ustar <- u[cbind(seq_len(nrow(V)), max.col(u, ties.method = "first"))] cnt <- cnt + ((ustar - rgumbel(nrow(V))) >= cut_true[KBUY - 1]) } ex_true <- cnt / 1000 m[r, ] <- c(sqrt(mean((hb_beta_i(f) - sp$est$Bmat_true)^2)), sqrt(mean((pred_exceed_panel(f, sp) - ex_true)^2)), hit_rate(f, sp), mean(f$draws$cut[, KBUY - 1]) - cut_true[KBUY - 1]) } res_all$nongumbel <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS)) knitr::kable(rbind(mean = res_all$nongumbel$mean, mc_se = res_all$nongumbel$se), digits = 4) ``` Preference recovery and demand prediction survive. The misspecification is absorbed almost entirely by `bias_c3`, the displacement of the top-2-box cut point. This is the pattern promised at the top: the cut points act as a sponge. ## 5. Correlated alternatives Appendix D proves the framework survives any GEV error structure once the inclusive value is redefined as $\ln G(e^V)$. Fitting the independent-Gumbel model to nested-logit data therefore misspecifies exactly one thing: it uses $\ln \sum_j e^{V_j}$ where the truth is $\ln G(e^V)$. Four alternatives in two nests at dissimilarity 0.5, which is a within-nest error correlation of 0.75. ```{r} #| label: nested NEST <- list(groups = c(1,1,2,2), lambda = c(0.5, 0.5)) m <- matrix(NA_real_, REPS, 4, dimnames = list(NULL, c("rmse_bi","rmse_ex","hit","bias_c3"))) for (r in seq_len(REPS)) { sp <- gen_panel(seed = 96000 + 10*r, nest = NEST) f <- fit_dual_hb(sp$est, mcmc = MCMC, seed = 97000 + r, verbose = FALSE) set.seed(98000 + r) d2 <- sp$holdout$design Bx <- sp$holdout$Bmat_true[d2$respondent_id, , drop = FALSE] V <- matrix(rowSums(d2$X * Bx), ncol = J, byrow = TRUE) cnt <- numeric(nrow(V)) for (mc in seq_len(1000)) { u <- V + rnested(nrow(V), J, NEST) ustar <- u[cbind(seq_len(nrow(V)), max.col(u, ties.method = "first"))] cnt <- cnt + ((ustar - rgumbel(nrow(V))) >= cut_true[KBUY - 1]) } ex_true <- cnt / 1000 m[r, ] <- c(sqrt(mean((hb_beta_i(f) - sp$est$Bmat_true)^2)), sqrt(mean((pred_exceed_panel(f, sp) - ex_true)^2)), hit_rate(f, sp), mean(f$draws$cut[, KBUY - 1]) - cut_true[KBUY - 1]) } res_all$nested <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS), nest = NEST) knitr::kable(rbind(mean = res_all$nested$mean, mc_se = res_all$nested$se), digits = 4) ``` Again the cut point absorbs it, and in the direction theory predicts. Nesting makes $G(e^V) < \sum_j e^{V_j}$, so the fitted model overstates the inclusive value, and the top-2-box cut compensates by moving upward. The sign of `bias_c3` is therefore a prediction, not a fitted artifact. ## 6. Task-order drift Respondents are not static across a survey. The outside good may drift as they work through it, which is the phenomenon documented in the no-choice literature. Two variants: a deterministic linear drift, and Hung-style probabilistic updating toward the previous task's inclusive value. ```{r} #| label: drift variants <- list(linear = list(type = "linear", delta = 0.06), update = list(type = "update", rho = -2.5)) drift_out <- list() for (vn in names(variants)) { m <- matrix(NA_real_, REPS, 4, dimnames = list(NULL, c("rmse_bi","bias_c3","order_slope","order_t"))) for (r in seq_len(REPS)) { sp <- gen_panel(seed = 95000 + 10*r, drift = variants[[vn]]) f <- fit_dual_hb(sp$est, mcmc = MCMC, seed = 96000 + r, verbose = FALSE) d <- sp$est$design mubar <- task_mubar(d, hb_beta_i(f)) chat <- colMeans(f$draws$cut)[KBUY - 1] resid <- as.numeric(sp$est$y >= KBUY) - (1 - plogis(chat - mubar)) tt <- rep(seq_len(T_EST), times = N) reg <- summary(lm(resid ~ tt)) m[r, ] <- c(sqrt(mean((hb_beta_i(f) - sp$est$Bmat_true)^2)), chat - cut_true[KBUY - 1], reg$coefficients[2,1], reg$coefficients[2,3]) } drift_out[[vn]] <- list(metrics = m, mean = colMeans(m), se = apply(m, 2, sd)/sqrt(REPS)) } res_all$drift <- drift_out knitr::kable(do.call(rbind, lapply(names(drift_out), function(v) data.frame(variant = v, t(round(drift_out[[v]]$mean, 4))))), row.names = FALSE) ``` The static model's cut points absorb the average drift while part-worths stay accurate. More useful for practice: a one-line diagnostic detects it. Regress the residual exceedance on task position, and under linear drift the $t$ statistic is large and consistently signed, so an analyst can find this in real data without fitting a dynamic model. ## 7. What the five have in common ```{r} #| label: summary knitr::kable(data.frame( misspecification = c("wrong link","heterogeneous cuts","normal errors", "correlated alternatives","task-order drift"), `preference recovery` = c("unaffected","unaffected","unaffected", "slightly worse","unaffected"), `demand prediction` = c("second-order","unaffected","unaffected", "unaffected","unaffected"), `where the damage lands` = c("tails only","nowhere material","cut points", "cut points","cut points"), check.names = FALSE), row.names = FALSE) ``` Four of the five leave the quantities practitioners actually use essentially intact, and deposit their damage in the cut points. The reason is structural rather than lucky: the forced choice identifies preferences and is untouched by any of these, while the cut points are the model's only free location and so absorb whatever the second stage cannot express. The implication cuts both ways. Demand estimates are robust. The *structural* reading of the cut points, as probability thresholds that give scale labels quantitative meaning, is the fragile part, and it is exactly what our model's probability-scale interpretation asks a reader to believe. ## 8. Results written ```{r} #| label: save PROJ <- if (file.exists("_quarto.yml")) "." else ".." OUT <- file.path(PROJ, "R", "output") dir.create(OUT, showWarnings = FALSE, recursive = TRUE) settings <- list(N = N, T_EST = T_EST, T_HOLD = T_HOLD, J = J, W = W, beta_bar = beta_bar, Sigma = Sigma_true, cut = cut_true, REPS = REPS, MCMC = MCMC, KBUY = KBUY) for (cell in c("wronglink","hetcut","nongumbel","nested","drift")) { saveRDS(list(cell = cell, res = res_all[[cell]], settings = settings), file.path(OUT, paste0("hb_simstudy_", cell, ".rds"))) } cat("wrote five hb_simstudy_*.rds files\n") ``` ## 9. Related material | where | what | |---|---| | notebook 05 | the tests that catch two of these before they do damage | | notebook 07 | why link discrimination depends on design, not sample size | | Appendix D | the GEV result behind the correlated-alternatives cell |