Skip to main content

Aggregate Recovery

Both links, three configurations, and forty independent datasets

Authors

Dan Yavorsky

Geoffery Zheng

Published

October 24, 2026

\[ \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

Notebook 01 worked one configuration of the model end to end. This one asks whether that was luck.

Three questions, in increasing order of what they would catch:

  1. Does recovery hold across configurations and both links? One cell showing recovery is an anecdote; six is a pattern.
  2. Is the estimator unbiased, and are its standard errors honest? Recovery on a single dataset says nothing about repeated sampling. Forty independent datasets per link do.
  3. Is the residual bias the ordinary finite-sample kind? One parameter comes in near the conventional threshold. Ordinary maximum-likelihood bias shrinks like \(O(1/n)\); a structural problem does not.

The machinery below is the general form of what notebook 01 built: it carries both links rather than one link alone, and it is defined once here rather than twice.

This notebook is the source of aggregate_recovery.rds, replication_study.rds and bias_scaling_check.rds.

In [1]:
set.seed(1)
options(digits = 5)

1. The machinery

The general form of what notebook 01 built. Two differences: these carry both links rather than one link alone, and they are defined once rather than twice.

In [2]:
Designs, the scale maps, and the cut-point parameterization 

# Standard Gumbel quantile function and CDF. Exact inverses; used in both
# directions. rgumbel draws by inverse transform.
alpha_to_cut <- function(a) -log(-log(a))   # probability -> utility,  F^{-1}
cut_to_alpha <- function(c) exp(-exp(-c))   # utility -> probability,  F
rgumbel      <- function(n) alpha_to_cut(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
}

# Unconstrained working parameters: (beta, c_1, log gaps). See notebook 01 for
# the derivation and the Jacobian used for delta-method standard errors.
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) c(cut[1], log(diff(cut)))

3. Repeated sampling: bias and standard-error calibration

One dataset per cell shows recovery. Forty independent datasets per link show whether the estimator’s sampling behavior matches the asymptotics: is it unbiased, and are the model-based standard errors the right size?

In [10]:
R_REPS <- 40
N_REP  <- 5000

rep_pars <- list(
  beta   = c(a2 = 0.8, a3 = -0.5, b2 = 0.4, b3 = 1.0, price = -0.9),
  cutB   = c(-1.0, 0.2, 1.2, 2.2),
  alphaA = c(0.10, 0.30, 0.60, 0.85)
)

run_model <- function(link) {
  cut_true <- if (link == "cloglog") alpha_to_cut(rep_pars$alphaA) else rep_pars$cutB
  truth <- c(rep_pars$beta, cut_true)
  est <- se <- matrix(NA_real_, R_REPS, length(truth))
  conv <- integer(R_REPS)
  for (r in seq_len(R_REPS)) {
    design <- make_design(N_REP, J, seed = 10000 + r)
    dat <- simulate_dual(design, rep_pars$beta, cut_true, link = link,
                         seed = 20000 + r)
    fit <- fit_dual_mle(dat)
    est[r, ] <- c(fit$beta, fit$cut)
    se[r, ]  <- c(fit$se_beta, fit$se_cut)
    conv[r]  <- fit$convergence
  }
  bias   <- colMeans(est) - truth
  emp_sd <- apply(est, 2, sd)
  data.frame(
    link = link,
    param = c(names(rep_pars$beta), paste0("c", 1:(W - 1))),
    truth = truth, mean_est = colMeans(est), bias = bias,
    t_bias = bias / (emp_sd / sqrt(R_REPS)), emp_sd = emp_sd,
    mean_se = colMeans(se), se_ratio = colMeans(se) / emp_sd,
    n_nonconv = sum(conv != 0))
}

rep_res <- do.call(rbind, lapply(c("cloglog", "logit"), run_model))
rownames(rep_res) <- NULL
knitr::kable(rep_res, row.names = FALSE, digits = 4)
link param truth mean_est bias t_bias emp_sd mean_se se_ratio n_nonconv
cloglog a2 0.8000 0.7922 -0.0078 -1.4765 0.0333 0.0338 1.0159 0
cloglog a3 -0.5000 -0.4986 0.0014 0.1908 0.0458 0.0410 0.8949 0
cloglog b2 0.4000 0.3872 -0.0128 -2.0506 0.0395 0.0392 0.9932 0
cloglog b3 1.0000 0.9925 -0.0075 -1.1189 0.0421 0.0370 0.8774 0
cloglog price -0.9000 -0.8971 0.0029 0.7299 0.0249 0.0266 1.0691 0
cloglog c1 -0.8340 -0.8364 -0.0024 -0.2915 0.0515 0.0578 1.1211 0
cloglog c2 -0.1856 -0.1969 -0.0113 -1.6206 0.0440 0.0520 1.1824 0
cloglog c3 0.6717 0.6680 -0.0037 -0.5130 0.0461 0.0509 1.1022 0
cloglog c4 1.8170 1.8150 -0.0020 -0.2292 0.0554 0.0534 0.9633 0
logit a2 0.8000 0.7923 -0.0077 -1.3293 0.0364 0.0383 1.0518 0
logit a3 -0.5000 -0.4984 0.0016 0.2203 0.0462 0.0447 0.9668 0
logit b2 0.4000 0.3866 -0.0134 -1.9221 0.0442 0.0432 0.9781 0
logit b3 1.0000 0.9926 -0.0074 -1.0356 0.0454 0.0413 0.9102 0
logit price -0.9000 -0.8961 0.0039 0.8903 0.0276 0.0298 1.0815 0
logit c1 -1.0000 -1.0084 -0.0084 -0.7730 0.0687 0.0668 0.9729 0
logit c2 0.2000 0.1980 -0.0020 -0.2155 0.0588 0.0608 1.0338 0
logit c3 1.2000 1.1992 -0.0008 -0.0778 0.0659 0.0602 0.9128 0
logit c4 2.2000 2.1923 -0.0077 -0.7889 0.0621 0.0630 1.0151 0 
In [11]:
cat(sprintf("max |t_bias| = %.2f over %d tests (5%% critical ~ 2.0; Bonferroni ~ 3.0)\n",
            max(abs(rep_res$t_bias)), nrow(rep_res)))
max |t_bias| = 2.05 over 18 tests (5% critical ~ 2.0; Bonferroni ~ 3.0)
cat(sprintf("standard-error calibration: mean ratio = %.3f (1 = perfect)\n",
            mean(rep_res$se_ratio)))
standard-error calibration: mean ratio = 1.008 (1 = perfect)

The calibration ratio compares the average model-based standard error to the empirical standard deviation of the estimates across the forty datasets. At one, the standard errors are telling the truth.

4. Is the residual bias just \(O(1/n)\)?

One parameter in the table above carries a \(t\) statistic near the conventional threshold. Two explanations compete: ordinary finite-sample MLE bias, which shrinks like \(1/n\), or something structural, which would not shrink. The test is to quadruple \(n\) on fresh seeds and see whether the bias falls by about four.

In [12]:
N_BIG <- 20000
cut_cloglog <- alpha_to_cut(rep_pars$alphaA)
truth_cloglog <- c(rep_pars$beta, cut_cloglog)

est_big <- matrix(NA_real_, R_REPS, length(truth_cloglog))
for (r in seq_len(R_REPS)) {
  design <- make_design(N_BIG, J, seed = 30000 + r)
  dat <- simulate_dual(design, rep_pars$beta, cut_cloglog, link = "cloglog", seed = 40000 + r)
  fit <- fit_dual_mle(dat)
  est_big[r, ] <- c(fit$beta, fit$cut)
}

bias_big   <- colMeans(est_big) - truth_cloglog
emp_sd_big <- apply(est_big, 2, sd)
bias_res <- data.frame(
  param = c(names(rep_pars$beta), paste0("c", 1:(W - 1))),
  truth = truth_cloglog, mean_est = colMeans(est_big), bias = bias_big,
  t_bias = bias_big / (emp_sd_big / sqrt(R_REPS)), emp_sd = emp_sd_big)
rownames(bias_res) <- NULL
knitr::kable(bias_res, row.names = FALSE, digits = 4)
param truth mean_est bias t_bias emp_sd
a2 0.8000 0.7983 -0.0017 -0.6370 0.0169
a3 -0.5000 -0.5007 -0.0007 -0.2031 0.0215
b2 0.4000 0.3955 -0.0045 -1.4826 0.0190
b3 1.0000 1.0005 0.0005 0.2039 0.0163
price -0.9000 -0.8986 0.0014 0.7377 0.0121
c1 -0.8340 -0.8324 0.0016 0.3857 0.0268
c2 -0.1856 -0.1863 -0.0007 -0.1654 0.0255
c3 0.6717 0.6714 -0.0003 -0.0766 0.0235
c4 1.8170 1.8157 -0.0012 -0.3501 0.0225 
In [13]:
b2_small <- rep_res$bias[rep_res$link == "cloglog" & rep_res$param == "b2"]
b2_big   <- bias_res$bias[bias_res$param == "b2"]
cat(sprintf("b2 bias at n = 5,000:  %+.4f\n", b2_small))
b2 bias at n = 5,000:  -0.0128
cat(sprintf("b2 bias at n = 20,000: %+.4f\n", b2_big))
b2 bias at n = 20,000: -0.0045
cat(sprintf("ratio %.1fx (O(1/n) predicts ~4x); max |t_bias| now %.2f\n",
            b2_small / b2_big, max(abs(bias_res$t_bias))))
ratio 2.9x (O(1/n) predicts ~4x); max |t_bias| now 1.48

The bias shrinks in proportion to \(n\), which is what the finite-sample explanation predicts and what a structural problem would not do.

5. The identification condition, violated

@prp-ident says the design matrix augmented with a column of ones must have full column rank. The canonical way to break that is to give every alternative the same extra attribute, so we append a constant column and refit. The theory predicts a specific signature rather than an outright failure: the likelihood stays perfectly flat along the direction that trades the constant’s part-worth against an equal shift in every cut point, so the optimiser still converges, but where it stops depends on where it started.

In [14]:

des_bad <- make_design(4000, J, seed = 12345)
rc_ok   <- id_rank_check(des_bad$X)

des_bad$X <- cbind(des_bad$X, const = 1)   # the deliberate violation
des_bad$P <- ncol(des_bad$X)
rc_bad    <- id_rank_check(des_bad$X)

pars     <- param_sets$set1_moderate
cut_true <- alpha_to_cut(pars$alphaA)
dat_bad  <- simulate_dual(des_bad, c(pars$beta, const = 0), cut_true,
                          link = "cloglog", seed = 909)

# Three starting values differing only in the constant's coefficient. The
# singular Hessian makes the standard errors NaN, which is part of the point.
starts <- c(-1, 0, 1)
fits <- suppressWarnings(lapply(starts, function(s0)
  fit_dual_mle(dat_bad, start = c(rep(0, des_bad$P - 1), s0, cut_to_par(cut_true)))))

bconst <- sapply(fits, function(f) f$beta[des_bad$P])
ridge_tab <- data.frame(
  start_const = starts,
  nll         = sapply(fits, `[[`, "nll"),
  beta_const  = bconst,
  c1_raw      = sapply(fits, function(f) f$cut[1]),
  `c1_-_const` = sapply(fits, function(f) f$cut[1]) - bconst,
  min_eig     = sapply(fits, `[[`, "min_eig"),
  check.names = FALSE
)

# Everything the proof says is identified: the non-constant part-worths and
# the differences cut_w - beta_const.
ident <- sapply(fits, function(f)
  c(f$beta[1:(des_bad$P - 1)], f$cut - f$beta[des_bad$P]))
ridge_res <- list(
  rank_ok     = rc_ok,
  rank_bad    = rc_bad,
  table       = ridge_tab,
  nll_spread  = diff(range(ridge_tab$nll)),
  raw_spread  = diff(range(bconst)),
  ident_spread = max(apply(ident, 1, function(r) diff(range(r))))
)

knitr::kable(ridge_tab, digits = 6)
start_const nll beta_const c1_raw c1_-_const min_eig
-1 9738.5 -0.511471 -1.32256 -0.81109 0e+00
0 9738.5 -0.011471 -0.82256 -0.81109 0e+00
1 9738.5 0.488529 -0.32256 -0.81109 -1e-06 
In [15]:
cat("rank of [X 1] without the constant:", rc_ok$rank,  "of", rc_ok$required,
    "->", if (rc_ok$pass) "identified" else "not identified", "\n")
rank of [X 1] without the constant: 6 of 6 -> identified 
cat("rank of [X 1] with the constant:   ", rc_bad$rank, "of", rc_bad$required,
    "->", if (rc_bad$pass) "identified" else "not identified", "\n\n")
rank of [X 1] with the constant:    6 of 7 -> not identified 
cat("spread across starts, raw beta_const:", signif(ridge_res$raw_spread, 4), "\n")
spread across starts, raw beta_const: 1 
cat("spread across starts, identified:    ", signif(ridge_res$ident_spread, 4), "\n")
spread across starts, identified:     1.448e-07 
cat("spread in log-likelihood:            ", signif(ridge_res$nll_spread, 4), "\n")
spread in log-likelihood:             6.912e-11

Appending the constant drops the rank by one, exactly as @prp-ident requires. The three fits reach the same log-likelihood to nine decimal places while their raw \((\beta_{\text{const}}, \mathbf{c})\) differ by the full distance between the starting values: the optimiser is sliding along a flat ridge and stopping wherever it entered. What the proof says is identified agrees across all three fits to six decimal places, which is the practical content of the proposition: the data pin down \(c_w - \beta_{\text{const}}\) and the non-constant part-worths, and nothing pins down where along the ridge those sit.

6. Results written

In [16]:
PROJ <- if (file.exists("_quarto.yml")) "." else ".."
OUT  <- file.path(PROJ, "R", "output")
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)

saveRDS(list(recovery = results,
             settings = list(N_TASKS = N_TASKS, J = J, W = W,
                             param_sets = param_sets)),
        file.path(OUT, "aggregate_recovery.rds"))
saveRDS(rep_res,  file.path(OUT, "replication_study.rds"))
saveRDS(bias_res, file.path(OUT, "bias_scaling_check.rds"))
saveRDS(ridge_res, file.path(OUT, "identification_ridge.rds"))

cat("wrote aggregate_recovery.rds, replication_study.rds, bias_scaling_check.rds,",
    "identification_ridge.rds\n")
wrote aggregate_recovery.rds, replication_study.rds, bias_scaling_check.rds, identification_ridge.rds