Default change. run_mxlogit() and run_mnlogit() now default to
scale_vars = "bhhh", below; run_nestlogit() gains scale_vars (its last
argument) and keeps "none". A fit that uses the default reaches the same
estimator's maximum through other coordinates, so its estimates, iteration
counts and convergence codes change. With any scale_vars but "none" (so
now by default for the mixed and multinomial logits), settings in parameter
units in control (nloptr's xtol_abs; optim's parscale, and pgtol for
the mixed and nested logits, which optim fits by L-BFGS-B) and custom
optimizers act on the scaled coordinates. scale_vars = "none" restores the
0.2.x path, for a mixed logit together with the 0.2.x start values (next
section); for a saved fit, update(fit, scale_vars = "none") (and, for a
mixed logit, theta_init).
scale_vars chooses how the optimizer's coordinates are scaled; it never
changes the estimator or the data. The likelihood is evaluated at the
parameters, on the design as supplied, while the optimizer works on rescaled
parameters, theta = m * theta_t + c, with the gradient rescaled by the chain
rule. Start values and bounds are given in natural units and mapped;
estimates and standard errors are reported in natural units.
scale_vars = "bhhh" scales each parameter by the inverse square root of
the diagonal of the BHHH (outer product of gradients) matrix at the start
values, rounded to a power of two, from one extra gradient pass. L-BFGS
starts every iteration from a single multiple of the identity as its
curvature, and in a model with many alternatives the curvatures of the
constants (large and small alternatives) and of the Cholesky entries spread
over orders of magnitude that no column scale reaches. On simulated designs
(census-style, claims-style and small) and on the shipped mode_choice
data, with the default optimizer (nloptr L-BFGS) and a budget of 2,000
likelihood evaluations, "bhhh" reached the maximum (log-likelihood within
1e-6 of a Newton-polished reference, Newton decrement below 1e-6) in (base
runs from the 0.2.x start values, which the mixed logit's default no longer
uses (next section); reruns from starts perturbed by 1e-13 moved the counts
by up to a quarter, never their order):
"none" and "sd" had not
arrived after 2,000 (they stopped up to 90 log-likelihood units short, and
one of those runs reported convergence);"none"
and 998 for "sd";"none";mode_choice included: 11 to 30, against 19 to 108 for
"none" and 20 to 31 for "sd"; on a small mixed logit "none" and
"sd" both failed at once.Over the test suite's default fits, those capped by a test's evaluation
limit included, "bhhh" took 20% fewer iterations for the multinomial logit
(1,532 against 1,915 over 149 fits) and 6% fewer for the mixed logit (3,591
against 3,832 over 146, from the 0.2.x start). The scales fix the
optimizer's metric, not a start that is wrong in a column's units: with the
census-style design's columns rescaled by 10^3 or 10^-3, no choice reached
the mixed logit's maximum in 2,000 evaluations from the 0.2.x start
(Cholesky diagonal 0.5 in every column's units), which is why the mixed
logit's default start is now in the columns' units (next section); from the
new start "bhhh" reaches the maximum in 54 evaluations ("sd" and
"none" still do not). The multinomial logit, whose start does not depend
on units, reached its own maximum on that design with "bhhh" in 38
evaluations, against about 1,300 with "sd", while "none" reported
convergence 1,920 log-likelihood units short. Where the diagonal is unusable
at the start values (no finite, positive entry, as when the likelihood
overflows there), the parameters are left unscaled, with a message.
What moves with the default. Where the likelihood is well identified,
estimates move within the optimizer's tolerance, and on the designs above
(the rescaled one from a start in its units) fits that stopped short of the
maximum, or reported a convergence they had not reached, reached it.
Standard errors are computed from the equilibrated information (below); they
move with the estimates, and where the information is near singular (a
variance at its boundary) they can differ widely. Where the data barely
identify a direction, estimates can move along it: on a 60-person panel in
our reference battery, the constants moved together by 1.0, 0.0005 standard
errors, at log-likelihoods 1e-7 apart. On small, degenerate likelihoods the
two paths can stop at different points: among the test suite's fits,
"bhhh" reached a log-likelihood 0.96 higher in one fit whose coefficients
diverge (the data separate) and stopped short in two (by 0.09 and 0.02; an
alternative never chosen, a random coefficient's variance going to zero).
With optimizer = "optim", BFGS for the multinomial logit (its stopping
rule reltol, relative to the log-likelihood, about 1.5e-8 of its value)
and L-BFGS-B for the mixed and nested logits (factr, about 2.2e-9) can
stop before the accuracy above, scaled or not: on the claims-style
multinomial logit, optim reported convergence 0.44 log-likelihood units
short with "bhhh" (23 evaluations), and 0.52 short with "none" (778
evaluations, with maxit raised to 2,000). control = list(reltol = ) or
list(factr = ) tightens them; maxit (100 by default) caps the
iterations.
run_nestlogit() keeps "none" as its default. Its "bhhh" scales are
fixed at the start values, and when the nest parameters end far from their
start (0.5) they can mislead the optimizer: over the test suite's nested
logit fits, "bhhh" took 14% more iterations than "none" (21,190 against
18,650 over 102 fits; one fit 540 against 29, to the same maximum), although
on one simulated design with a covariate in large units "none" failed
where "bhhh" converged in 27 evaluations, and "sd" in 35.
"sd", "mad" and "iqr" keep their meaning (the optimizer works on the
parameters the design would have with its columns divided by their scales)
but take the same route: the data are never divided, and run_mnlogit() no
longer holds a scaled copy of X while it fits (on the claims-style
design's 8.7 million rows and 10 fixed-coefficient columns, a scaled fit's
peak heap fell from 1.39 to 0.10 GB, close to an unscaled fit's 0.08 GB).
Their estimates move within 1.7e-5 standard errors on our reference battery
(1.2e-6 over the test suite's fits).
The standard errors of a scaled fit are computed in natural units, at the
natural estimates, and its information matrix is equilibrated before it is
inverted: scaled by the powers of two nearest (on a log scale) the inverse
square roots of its diagonal's magnitudes, which adds no rounding. A fit
used to invert the information of the scaled design and transform the result
back, while vcov(fit, type = ) inverted the natural information as it was;
both lost digits when the design's columns span orders of magnitude (post
hoc, the inversion failed as computationally singular on such a design; for
one of the tests' mixed logit robust variances, whose information had
condition number 1.2e12, the standard errors computed raw were off by 2.1e-7
relative against an exact computation, equilibrated by 9e-16). A mixed
logit's variance at fit time is computed by the route vcov(fit, type = )
takes post hoc ("robust" for se_method = "sandwich"), with the fit's own
draws, so the two agree, bit for bit on one thread (draws supplied through
eta_draws are not regenerated post hoc). Unscaled fits now take the same
route (see "Every fit's information is equilibrated before it is inverted",
below).
The extra gradient pass of "bhhh" counts in the reported optimization
time. A scaled fit's theta_init must have one finite value per parameter;
a shorter or longer vector used to be recycled.
Fits record the map in param_scale and param_shift (1s and 0s for
"none"); sX and sW hold the column scales under "sd", "mad" and
"iqr", and 1s otherwise.
theta_init = NULL, run_mxlogit() now starts
each normal random coefficient at L_pp = 1 / h_p on the Cholesky diagonal,
where h_p is the typical step of its column: the distance between its values
for a column with two (a dummy: 1), the mean distance from a value that
fills more than half of the column over the other rows (a mostly-zero
column), and its standard deviation otherwise. It started at 0.5 whatever
the column's units. At the start, a step of h_p in a column now moves the
utility by a random amount with standard deviation 1 (the logit error's is
about 1.28), whether a price is in cents or in thousands. Rescaling a column
divides its row of the Cholesky factor (and its mean, with rc_mean = TRUE)
by the same factor and changes nothing else, and the start now moves with
it: in exact arithmetic, with scale_vars = "sd", "mad" or "iqr" the
optimizer's problem is the same in any units, and with "bhhh" the same up
to the power-of-two rounding of its scales. A log-normal coefficient keeps
L_pp = 0.5, the standard deviation of the logarithm of its random part,
which has no units; that random part has median one in the covariate's
units, though, so a log-normal covariate still belongs in units where its
coefficient is of order one. The estimator is unchanged: a fit that uses the
default start reaches the same maximum within the optimizer's tolerance,
except where one start or the other stops short of it, or the likelihood has
another local maximum or its supremum on a boundary (below). A theta_init
of zeros with log(0.5) on the Cholesky diagonal restores the 0.2.x start,
and with scale_vars = "none" the 0.2.x path.scale_vars, above) cannot repair a
start in the wrong units. On simulated designs whose columns we rescaled by
up to 10^±3 (census-style with 20,000 students, a cross-section of 10^6
rows) and on rescaled versions of the shipped mode_choice data, the
default "bhhh" from the 0.2.x start stopped 233 to 1,730 log-likelihood
units short of the maximum in 4 of 10 rescalings, three of them reporting
convergence (the fourth ran out of the budget), and needed 167 evaluations
in a fifth, against 42 on the original units. From the new start it reached
the maximum in all 10: in 0.84 to 1.09 times the evaluations it took on the
original units in 9 of them (54 on the census-style design rescaled by 10^3
and 10^-3, against 56 on its original units), and in 1.8 times as many in
the tenth, a rescaling of mode_choice (58 against 32; "bhhh"'s scales
round, and its stopping rule stops, differently in other units, while under
"sd" every rescaling of mode_choice took 29). Over these and the
original designs, 16 in all (claims-style with 518 parameters and
census-style with 200,000 students included; nloptr, 50 draws, a budget of
2,000 evaluations, the accuracy criterion above), the new start took at most
1.35 times the evaluations of the best of the five starts we compared on
each; the 0.2.x start missed 4 and took up to 3.9 times elsewhere. Against
the 0.2.x start: 157 against 219 on the claims-style design (213 in a
replicate; 198 in the study above), 44 against 49 on the census-style design
at 200,000 students, 56 against 42 at 20,000. Every column of these designs
has many values, so their step is the standard deviation."bhhh" fixes its scales at the start, and
with too little heterogeneity there they are too generous for the Cholesky
diagonal, whose curvature grows with L_pp: from 0.25 / sd(W_p) "bhhh"
reported convergence on the cross-section 554 log-likelihood units short.
"none" and "sd" have no such scales, but from 0.5 / sd(W_p) both failed
at once on that design, where from 1 / sd(W_p) both reached the maximum (in
43 evaluations)."bhhh", reached a higher log-likelihood (up
to 0.81 higher, on a 30-situation design where the 0.2.x start reported
convergence short of it); five, all with scale_vars = "none" on
deliberately badly scaled designs, a lower one (0.18 to 0.27 lower): there
the unscaled optimizer can drive a random coefficient's variance toward
zero, where the likelihood flattens and its gradient vanishes. On a weakly
identified design in our reference battery whose supremum lies on the
boundary (a degenerate covariance), "none" and "sd" now stop at the
interior maximum that "bhhh" reached before and still reaches, 0.021 below
the supremum: the second random coefficient's variance given the first
(L_22²) moves there from nearly zero to about 1.4.lower/upper, it is moved to the
nearest bound, with a message naming the parameters; with the default
optimizer such a call used to stop with an error (nloptr refuses a start
outside its box). A supplied theta_init is used as given. Bounds that hold
no start (lower above upper, lower = Inf, upper = -Inf) now stop
run_mxlogit() before the fit, naming the parameters, whatever the start;
before, nloptr stopped with its own error, while optim's L-BFGS-B returned a
fit with convergence code 52 for lower above upper and treated an
infinite bound on the wrong side as no bound.stats::sd() a column at a time would copy each
column and an index of its rows: 1.2 GB more at 98 million rows).run_mxlogit() / prepare_mxl_data() gain person_col: the frequentist
panel mixed logit of Revelt and Train (1998). When supplied, all choice
situations of a decision maker share one draw of the random coefficients,
and the likelihood integrates the joint probability of that person's
choices over the taste distribution, instead of integrating every choice
situation separately as before. id_col must still identify choice
situations uniquely across the whole data set (ids do not restart within
person). NULL (the default) keeps every choice situation its own decision
maker, reproducing the existing cross-sectional likelihood.
simulate_mxl_data() gains T (choice situations per decision maker;
default 1, a cross-section) and records the realized random coefficients
in true_params$gamma_i (K_w x N), the realized tastes that
conditional_tastes()'s conditional means can be compared with (they
track them with shrinkage). T = 1 reproduces earlier versions' data
exactly.conditional_tastes() (class choicer_tastes,
method for choicer_mxl): the mean and SD of a decision maker's random
coefficients conditional on the choices they made (Revelt and Train 2000;
Train 2009, ch. 11) — the frequentist analogue of the hierarchical Bayes
beta_i summaries of run_hmnlogit(). Its print() reports, per
coefficient, the share of population taste variance revealed by choices
and the residual share (law of total variance), which should sum to
about one — a numerical check on convergence and the draw count S, not
a specification test (both identities are first-order conditions of the
likelihood, so they hold at a converged fit whatever the true taste
distribution).person_col, scores are per decision maker, so
se_method = "sandwich" / vcov(type = "robust") is already robust to
within-person dependence; cluster_col / cluster labels must then nest
decision makers (be constant within each). Weights become decision-maker
weights, constant within person (the objective is sum_n w_n log L_n).
WESML (choice-based) weighting is not supported with person_col: a
choice-based weight varies with the alternative a person chose, which can
differ across their situations.n_persons (the number of
decision makers; printed as "Respondents"), data$Ti (choice situations
per decision maker, in prepared order), and data$person_ids. nobs()
is unchanged — it still counts choice situations, as for HMNL.inst/simulations/mxl_panel_simulation.R, in the
style of mxl_simulation.R: simulates a panel, fits it with
person_col, reports recovery_table() and conditional_tastes(), and
contrasts standard errors against a cluster-robust cross-sectional fit of
the same data.scale_vars (see
"Scaling the optimizer's coordinates"), default (cross-sectional)
behavior is unchanged up to floating-point reassociation (<= 1e-10
relative in the estimation kernels, measured on the reference battery of
test configurations); prediction kernels are bit-identical (see also
"Mixed logit — prediction at population scale").NULL empty Ti is an error, and a
primary-thread memory check now stops when a single decision maker
stacks so many alternative rows (tens of millions) that its design rows
would need more than 2 GiB of scratch in every thread — a signal that
person_col identifies markets or some other high-cardinality grouping
rather than decision makers.weights vector and an
unnamed post-hoc cluster vector are rejected (with a pointer to
weights_col / id-named labels) whenever that order differs from id
order, and predict() / logsum() / consumer_surplus() note when
newdata lacks the person column (rows are then ordered by id rather than
in the in-sample, decision-maker-first order).conditional_tastes() on other fits errors informatively, pointing
hierarchical Bayes users to their beta_i summaries.person_col to use the within-person variation the
panel provides.S. The draw loop also sweeps a decision maker's rows in
memory order and no longer allocates per decision maker. On a synthetic
claims panel with heavy-tailed histories (up to 500 choice situations per
decision maker), peak working memory at S = 1000 fell from 6.4 GB to
0.12 GB and gradient evaluations ran 40% faster; on a 10^7-row panel they
ran 12-25% faster at S = 100-500. The analytical Hessian no longer
holds a decision maker's rows-by-draws matrix either: its peak working
memory on the claims panel at S = 500 fell from 1.7 GB to 0.2 GB. A
decision maker is split into batches only when its stacked alternative
rows times S exceed 2^18 (a long panel, a choice set of hundreds of
alternatives, or a very large S).X beta + W mu + delta) inside the parallel loop instead
of for all rows up front, so a likelihood evaluation no longer allocates
anything the length of the stacked design (what remains is 8 bytes of row
offsets per choice situation, and in a panel 8 per decision maker). On
synthetic 10^8-row panels, peak working memory fell from about 3 GB to
under 0.2 GB and an evaluation ran 19% faster at S = 25 (7-9% at
S = 100); fits at S = 50 on 10^7 rows took 9% less time per
evaluation. The kernels' offsets into the stacked design are now 64-bit.draws = "generate") are now formed in two passes
over each block: all the Halton uniforms first, then the inverse normal
CDF. The first (base-2) dimension is stepped by an exact integer odometer,
and the digits of the next ten (bases 3 to 31) are taken by division by a
compile-time constant. The draws are bit-identical to before on 64-bit
platforms, so this change by itself alters no estimate, standard error,
prediction or other post-estimation quantity; it lowers the cost of
regenerating the draws at each likelihood evaluation. With three random
coefficients and S = 100, gradient evaluations ran 32-36% faster on a
synthetic four-alternative cross-section (10^6 rows at 1 and 11 threads,
10^7 rows at 11), where regenerating the draws had taken about half the
time; 10-14% faster on a school-census-like design (thousands of schools,
10-100 per choice set); and showed no measurable change (under 3%, within
run-to-run noise) on claims-like panels with 20-200 hospitals per choice
set, whose cost lies elsewhere. On the four-alternative data a
generate-mode evaluation now takes about 1.35 times as long as a
store-mode one, down from about 2 times. The generator's
digit-permutation tables are also smaller: 1.3 MB instead of 11 MB at the
128-dimension maximum.scale_vars (see "Scaling the
optimizer's coordinates").run_mnlogit() and run_nestlogit()
(likelihood, gradient, Hessian, BHHH and scores) and behind their
post-estimation methods (predictions, shares, elasticities, diversion
ratios, blp_contraction() and nl_blp_contraction()), whose
contractions also form X beta once instead of at every iteration. The
nested logit kernels also compute each choice situation's probabilities
in buffers allocated once per thread, where they built about a dozen
Armadillo vectors per situation. On a synthetic design of 17.5 million
rows (500,000 choice situations of 20 to 50 alternatives, 10 covariates),
a multinomial logit gradient evaluation took 23% less time at one thread
and 39% less at eleven, moved 0.5 GB less through the heap, and peaked at
0.14 GB of working memory instead of 0.57 GB; twenty iterations of the
BLP contraction took 42% and 70% less time (medians of six alternated
runs). On 7 million rows (200,000 choice situations of 20 to 50
alternatives in 10 nests), a nested logit gradient evaluation took 37%
less time at one thread and 46% less at eleven, and peaked at 0.06 GB of
working memory instead of 0.23 GB; twenty iterations of its
contraction took 47% and 65% less time and moved 0.1 GB through the heap
instead of 10.9 GB (medians of three alternated runs; with 40 nests, 42%
and 50% for the gradient and 51% and 64% for the contraction).
Estimates, standard errors, predictions and every post-estimation
quantity are unchanged, bit for bit.input_data list or a direct call can supply, in a defined way. A missing
or negative alternative code is reported as such, where it went through an
undefined conversion (and read as 0 on ARM processors), and a missing nest
code always as an invalid one; a missing or negative choice is an error,
where on ARM, with an outside option, it counted as choosing that option;
and an index beyond the integer range becomes NA, with R's warning, and is
then treated as a missing one. blp_contraction() and
nl_blp_contraction() given no choice situations and an empty delta now
stop at their check of the weights, where they read past the empty vector
or, with an outside option, failed an Armadillo bounds check.predict() (probabilities and shares),
logsum(), elasticities(), diversion_ratios(), blp() and the exported
mxl_blp_contraction() (and so consumer_surplus(), gof() and the fit
statistics of summary()) now read the stacked design's alternative codes
in place and simulate each choice situation in buffers allocated once per
thread. They used to copy the codes twice on every call and allocate about
ten objects per choice situation, among them its rows-by-draws matrix of
random utilities and, for the log-sums, a temporary at every draw; with
stored draws, the threads also took turns at a lock shared by all of them
for every situation's draws. The BLP contraction also forms X beta and
W mu once instead of at every iteration. The working arrays, among them
the alternative-by-alternative accumulators of elasticities and diversion
ratios, are allocated before the parallel loop, so running out of memory now
stops with an R error that says how much each thread needs and suggests
fewer threads, instead of possibly aborting R from inside the loop. Results
are unchanged, bit for bit, with a BLAS whose results do not depend on where
its operands sit in memory, such as the reference BLAS or OpenBLAS. That
holds for each choice situation's probabilities, mean utilities and log-sum
at any number of threads, and for sums over situations (shares,
elasticities, diversion ratios, and so BLP inversions) at one thread; at
more threads such sums vary in their last bits with the order in which the
threads take the situations, as before. With other BLAS libraries, results
may differ from earlier versions' in their last bits.
S = 100), at one and eleven threads, with draws generated
on the fly and, in parentheses, stored (medians of three alternated runs).
Probabilities took 16% and 20% less time (23% and 39%) and moved 0.16 GB
through the heap in 51 allocations instead of 11.8 GB in 2.2 million.
Log-sums took 27% and 26% less time (65% and 90%) and moved 0.08 GB in 19
allocations instead of 18.8 GB in 22 million. Shares took 14% and 16% less
time (20% and 36%) and moved 0.08 GB in 49 allocations instead of 11.4 GB
in 2.0 million. Five iterations of the BLP contraction took 14% and 17%
less time (17% and 29%) and moved 0.08 GB instead of 67.5 GB. On a
claims-like panel of similar size, the same kernels took 14-27% less time
for the probabilities, 24-44% less for the log-sums and 12-23% less for
the shares and the contraction.NA, with R's warning, and is then treated as a missing one.
With draws generated on the fly, an alternative-level W with fewer rows
than the largest alternative code is now reported as with stored draws,
where it stopped at an Armadillo bounds error.draws = "store", the default), predict(), logsum(),
consumer_surplus(), elasticities(), diversion_ratios(), blp() and
gof() (and so the fit statistics of summary()) integrate each choice
situation over its own block of Halton draws, and used to build the blocks
of all N situations as one K_w x S x N cube on every call (for a panel fit
too: predictions are unconditional, one block per situation). They now form
each situation's block on the fly, in the thread that simulates the
situation, and hold no cube. choicer's own Halton generator, without its
digit permutations, computes the points of randtoolbox::halton() digit by
digit as randtoolbox does, and R's qnorm() maps them to normal draws as
get_halton_normals() does. The draws, and so the results, are the cube's
bit for bit where randtoolbox and choicer are compiled with the same
floating-point contraction, as on CRAN's platforms: per choice situation at
any number of threads, and sums over situations at one thread. Otherwise
some of the uniforms differ in their last bit, which moves those draws by at
most a few parts in 10^9, in the far tails. Together with the kernel changes
above, on the census design at 98 million rows (2.45 million choice
situations, whose cube takes 5.9 GB), predict() of choice probabilities
with stored draws took 9.0 s instead of 76 s, moved 1.6 GB through the heap
in 2,431 allocations instead of 151 GB in 24 million, and its working memory
peaked at 1.6 GB, the probabilities and mean utilities it returns, instead
of 14.4 GB; on the claims-like panel at 85 million rows (1.5 million choice
situations, whose cube takes 3.6 GB), 6.4 s instead of 24 s and 1.4 GB
instead of 12.2 GB (medians of three alternated runs at eleven threads,
under background load). Without the cube these functions no longer stop past
store mode's 2^31 - 1 points (S times the number of choice situations),
which a panel fit could reach here although its estimation draws, one block
per decision maker, stayed within them. With more than 128 random
coefficients (the generator's table of primes) they build the cube, as
before.blp() on a store-mode fit used to build the cube once per call and read it
at every iteration of the contraction; it now forms the draws on the fly at
every iteration, in parallel, as predict() does. At eleven threads, each
evaluation of the shares (one per iteration, and one before the first) took
about 8.5 s on the census design and 6.2 s on the claims-like panel, against
8.3 s and 7.1 s before, when blp() also spent about 60 s and 20 s building
the cube once per call (slopes between one and seven iterations, single runs
under background load); its working memory peaked at 0.8 GB instead of
12.1 GB, and at 0.7 GB instead of 10.9 GB. The exported
mxl_blp_contraction() accepts gen_scramble = 2 for these store-mode
points formed on the fly, and now rejects other codes than 0, 1 and 2,
which it treated as 0.blp() gains keep_draws_bytes, the most memory it may use to keep the
draws across the iterations of its contraction (default 2^32 bytes,
4 GiB). When the draws of all N choice situations take at most that
(8 * K_w * S * N bytes), the first evaluation of the shares forms them, in
parallel, and the later ones reuse them; above it, every iteration forms
them again, as before. 0 never keeps them, and Inf keeps them whenever
they are formed on the fly and can be allocated; the deltas are the same
either way. The memory is taken per call, and where it is committed only
when used (Linux, macOS, containers and batch jobs with memory limits), a
budget above the memory actually free can end the R session rather than fall
back. At eleven threads, keeping them took each evaluation of the shares
from 8.5 s to 7.2 s on the census design with stored draws and from 9.4 s to
6.9 s with generated draws (9.1 s before this release; its 5.9 GB of draws
are above the default, so these runs used keep_draws_bytes = Inf), and on
the claims-like panel, whose 3.6 GB are within the default, from 6.4 s to
5.2 s with stored draws and from 7.0 s to 5.6 s with generated draws (8.4 s
before); single runs under background load, two alternated ones for the
claims-like panel's stored draws. Working memory then holds the draws: it
peaked at 6.7 GB and 4.3 GB.run_mxlogit(draws = "store") (the convenience workflow) now warns, before
it builds its draws, when they would take more than 1 GiB (8 * K_w * S * U
bytes, for the U decision makers of a panel fit or the choice situations of
a cross-section): the cube is held while the model is fitted, and
vcov(type = ), wesml_vcov() and conditional_tastes() rebuild it whole.
draws = "generate" stores no draws; with scramble = "none" it uses the
same Halton points, so its likelihood is store mode's up to rounding error
and its estimates agree to the optimizer's tolerance. The warning does not
change the fit; under options(warn = 2), which turns warnings into errors,
such fits now stop before estimating.run_mxlogit()'s
default standard errors, vcov(type = "hessian"), and the bread of the
robust, WESML and clustered variances) now works, for each decision maker,
only on the parameters that decision maker's likelihood involves: the
fixed coefficients, the means and Cholesky factor of the random
coefficients, and the constants of the alternatives that decision maker
faced; each choice situation's terms involve its own alternatives'
constants only. It used to carry every alternative's constant through
every decision maker, so that with thousands of alternatives most of its
work multiplied zeros; the gain is largest where each decision maker faces
a small share of the alternatives. On a synthetic school-census design with
2,620 parameters (2,602 schools, 10 to 100 per choice set, most students
with one choice set and some with up to three, three correlated random
coefficients, S = 100), where a student's likelihood involves 59 of the
2,620 parameters on average, the Hessian for 2,000 students took 0.35 s
instead of 24.8 s at eleven threads (1.9 s instead of 101 s at one), and
for 20,000 students 3.2 s instead of 228 s; its peak working memory fell
from 1.9 GB to 0.4 GB. On a claims-like panel with 500 hospitals (518
parameters, 20 to 200 hospitals per choice set, up to 500 visits per
patient) it took 3.4 s instead of 7.1 s for 1,500 patients and 26 s
instead of 55 s for 15,000 (medians of three alternated runs under
background load): a patient's visits involve about 60 of the hospitals on
average, and the work on each visit's alternatives, which this change
leaves as it was, now dominates.se_method = "bhhh", vcov(type = "bhhh") and, with squared weights, the
meat of the robust and WESML variances) now adds each decision maker's
rank-one update over the same parameters; it used to form every decision
maker's update over all of them. On the school-census design the BHHH
matrix for 20,000 students took 0.15 s instead of 25.1 s at eleven
threads, and for 200,000 students 1.4 s instead of 246 s; on the
claims-like panel, 15,000 patients took 0.88 s instead of 1.14 s (medians
of three alternated runs).vcov(type = "robust"), vcov(type = "cluster"), wesml_vcov()) and the clustered
variance of run_mxlogit(se_method = "cluster") no longer form the matrix
of every decision maker's score, one row per decision maker and one column
per parameter: 42 GB for 2 million students with 2,620 parameters, held
twice while it was weighted. (se_method = "sandwich" at fit time already
took its meat from the BHHH matrix with squared weights.) The robust meat
is now that BHHH matrix, the clustered meat sums the weighted scores
within clusters in C++ (a new internal kernel that splits a large
cluster's sum across the threads), and vcov(type = "bhhh") uses the BHHH
kernel as the fit does. On the school-census design with 200,000 students
at eleven threads, the meat of the robust variance took 1.4 s and 0.4 GB
of working memory instead of 12.2 s and 8.4 GB; the clustered meat 1.6 s
instead of 4.6 s with 2,000 clusters and 1.1 s instead of 3.2 s with five
(medians of three alternated runs). These variances also build their
store-mode draws once per call instead of twice. The cluster labels they
accept, their checks and messages, and the absence of a small-sample
correction are unchanged. The internal score kernel now refuses to form a
matrix above 2 GiB.options(choicer.max_memory = <bytes>) sets the limit instead, such as a
container's or a cluster job's, and Inf lifts the check
(?set_num_threads). Results are unchanged.se_method = "bhhh" is unchanged bit for bit (vcov(type = "bhhh") after the fit, now computed from it, by at most 2e-16), the
Hessian differs by at most 3.5e-15 relative (not at all on designs with
many alternatives), and the robust and clustered meats by at most 7e-16.
As before, at more than one thread, sums over decision makers vary in
their last bits with the order in which the threads take them; a decision
maker whose utilities overflow is left out of the Hessian and spreads
NaN over its parameters' rows and columns in the BHHH matrix and the
meats; and a non-finite term of a decision maker's Hessian (an
overflowing utility at a draw of zero weight) makes the rows and columns
of the parameters concerned NaN.get_halton_normals() builds the draw cube of draws = "store", at fit
time and again for post-estimation. It now fills the cube a block of units
at a time; before, it generated the whole Halton sequence in one call and
copied it into the cube unit by unit in an R loop. The draws are
bit-identical wherever the old code was well defined. For 2 million draw
units at S = 100 with three random coefficients (a 4.8 GB cube),
building the cube took 21 s instead of 27 s, peak memory fell from 3.0 to
1.4 times the size of the cube (the rest is R's garbage-collection slack),
and 1.3 thousand allocations replaced 8 million. The single call also
formed positions in its output with C int arithmetic, which is undefined
past 2^31 - 1 values (a 17 GB cube); each block now stays far below that.
The function now checks that S, N and K_w are positive whole
numbers, and stops before allocating a cube that store mode cannot
handle: more than 2^31 - 1 points (S * N), the largest starting index
randtoolbox::halton() accepts, where the old code failed with an
unrelated error. draws = "generate" has no such limit. choicer now
requires randtoolbox 1.31.0 or later, the release that added
halton(start = ).prepare_mnl_data(), prepare_mxl_data() and prepare_nl_data() (and so
run_mnlogit(), run_mxlogit() and run_nestlogit()) copy only the
columns the model uses from a data.frame, tibble or data.table, where they
copied the whole data set; they scan those columns for missing values one
at a time instead of through a rows-by-columns is.na() matrix, and check
that the covariates are numeric without copying them. prepare_nl_data()
no longer makes a second full copy to read the alternative-to-nest map. On
a synthetic claims-style panel of 8.7 million rows, peak heap use fell from
5.7 to 4.3 times the size of the input for prepare_nl_data() and from
4.5 to 4.0 times for prepare_mxl_data(); the saving grows with the
number of columns the model does not use. Apart from the storage change
below, the prepared objects are unchanged on our reference battery and on
every call the test suite makes.X and W in
C++ straight from the caller's covariate columns once the rows are
filtered and sorted. X and W are now always double precision: a
design matrix whose columns were all integer used to be an integer
matrix, which the estimation kernels converted to double on every call.
Estimates, standard errors and predictions are unchanged, bit for bit.qr(X, tol = 1e-7) moves past its rank) now factors a design of a
million or more values by row chunks and applies qr()'s own rank rule
to the resulting factor, which has the same cross-products. It drops the
same columns on our reference battery, on every call the test suite makes
and on a sweep of near-collinear designs; the two can disagree only on a
column whose dependence lies within rounding of the tolerance, where
qr()'s own answer already changes with the order of the rows. The check
no longer makes qr()'s copies of the design (about three times its
size), nor stops at its limit of 2^31 - 1 values, which ended the
preparation of any larger design (for instance 2.2 * 10^8 rows and 10
covariates) with "too large a matrix for LINPACK". Such designs can now
be prepared for every model.prepare_mxl_data() moved
through the heap on the 8.7-million-row claims-style panel; the
per-situation value is read from each situation's first row. With all the
changes above, on that panel and on a census-style panel of 9.8 million
rows, prepare_mxl_data() now takes 2.9-3.7 s instead of 5.7-7.6 s,
allocates 7.5-8.8 GB instead of 97-149 GB, and its peak heap use is 1.9
times the size of its input instead of 4.0-4.5 times; for
prepare_mnl_data() and prepare_nl_data() it is 2.0-2.1 times instead
of 4.4-5.7 times. Peak resident memory, which also counts memory R has freed
but not returned, is 2.3-2.4 times the input for all three instead of
5.3-8.1 times (medians of three runs).run_mxlogit(scale_vars = ) no longer holds a scaled copy of the design. It
used to divide X and W by their column scales with sweep() and keep
the results beside the original matrices for the whole fit, which also took
two design-sized temporaries to form. The scales themselves came from
apply(), which first copies the whole matrix. Scaling is now a change of
the optimizer's coordinates (see "Scaling the optimizer's coordinates"
above), so the data are never divided, and the scales are computed a column
at a time, collecting each column's temporaries before the next on large
designs. On the claims- and census-style panels above (8.7 and 9.8 million
rows, 13 covariates), the memory a scaled fit holds while it optimizes fell
by 0.9 and 1.0 GB, one copy of X and W. Fitted from the long data,
preparation included, a scaled fit's peak heap is now an unscaled one's: it
fell from 3.3 to 1.9 times the size of the long data (from 3.4 and 3.9 GB to
2.0 and 2.2 GB), with scale_vars = "sd" and "mad" alike. With the data
prepared beforehand, "sd" raises the peak by 0.06 GB at most (it raised it
by 1.9 and 1.6 GB) and "mad" by 0.3 GB at most (2.6 and 2.3 GB before).
Computing the scales adds 0.5-0.6 s to such a fit with "sd" and 2.8-3.3 s
with "mad" (1.3-1.5 s and 3.7-4.2 s before, with the scaled copies), and
"bhhh"'s extra gradient pass 0.35-0.46 s (single runs of one likelihood
evaluation and the BHHH variance, under background load; the earlier times
were medians of three alternated runs). The collection
after each column takes longer the more distinct character strings a session
holds, the data's own included: about 0.05 s with 10^7.
run_mnlogit(scale_vars = ) no longer holds a scaled copy of X either
(see above).prepare_mnp_data(), prepare_hmnl_data() and prepare_hmnp_data() (and
so run_mnprobit(), run_hmnlogit() and run_hmnprobit()) now also copy
only the columns the model uses, scan them for missing values one at a
time, and check that the covariates are numeric without copying them. The
hierarchical preparations also scan for non-finite covariate values a
column at a time, and look up the choice situations to drop only when a
row is flagged, where they grouped every row by choice situation twice.
Apart from the changes below, the prepared objects are unchanged on every
call the test suite makes and on a battery of edge cases; a used column
that shared its name with a working column of the old code (HAS_NA,
TASK_HAS_NA, HAS_BAD, TASK_HAS_BAD) stopped the preparation with an
error, and now works like any other.X in C++ from the caller's covariate columns,
as the MNL, MXL and NL preparations do; prepare_mnp_data() writes each
choice situation's differences from its base alternative in the same
pass. X is now always double precision there too (posterior draws are
unchanged, bit for bit). The hierarchical preparations also read the
alternative-level covariates by source row, call a covariate constant
within choice situations when every row equals its situation's first row,
and look for an alternative listed twice in one pass over the keys; both
checks used to evaluate R code for every choice situation. A few inputs
now prepare differently. The choice-situation id used as a covariate (a
task-order term, say) is now recognized as constant within every
situation: it is kept with a message, or dropped with a warning when there
is no outside option, where it used to be kept silently. Inputs whose
column names clashed with the preparations' own names are covered under
Corrections. Integer covariates whose differences or within-situation
ranges overflow 32-bit integers no longer fail. On synthetic inputs of
10 million rows in which the model uses every column, prepare_mnp_data()
(2.5 million choice situations of four alternatives) now takes 5.0 s instead
of 13.8 s, allocates 6.4 GB instead of 9.0 GB, and its peak heap use is
3.0 times the size of its input instead of 4.9 times; prepare_hmnl_data()
(500,000 respondents with five choice situations of four alternatives
each) takes 14.5 s instead of 31.7 s, allocates 12.2 GB instead of
15.5 GB, and peaks at 3.1 times its input instead of 4.9 times. Peak
resident memory is 4.5 times the input instead of 8.0 times for
prepare_mnp_data() and 4.7 times instead of 7.6 times for
prepare_hmnl_data() (medians of three runs).data.table::rowidv(); prepare_mnp_data() has none to number);
they used to evaluate R code for every choice situation. On the inputs of
the previous item this halves the time prepare_mnp_data() takes and
cuts that of prepare_hmnl_data() by two fifths; on the claims- and
census-style panels above, prepare_mnl_data(), prepare_mxl_data() and
prepare_nl_data() take 9-15% less time. Memory traffic grows by up to 3%
(data.table's working arrays for the sums), and peak heap use is unchanged
or lower (medians of three alternated runs). The prepared objects are
unchanged on every call the test suite makes and on a battery of edge
cases.optimizer = "optim", run_mnlogit(), run_mxlogit() and
run_nestlogit() no longer repeat the likelihood pass for a gradient at the
point just evaluated. stats::optim() asks for the objective and its
gradient through separate functions, and each call ran the pass, which
computes both, so a point at which optim asked for both was evaluated twice;
the gradient now comes from the pass that computed the objective. L-BFGS-B,
which fits the mixed and nested logits, asks for both at every point, so
those fits take about half the evaluations; BFGS, which fits the multinomial
logit, also asks for the objective alone, at its line search's trial points,
and saves one evaluation per gradient. Over the test suite's fits refitted
with optimizer = "optim", the evaluations (the "bhhh" pass included)
fell from 14,104 to 7,087 for the mixed logit (198 fits), from 12,941 to
6,470 for the nested logit (95) and from 4,671 to 3,118 for the multinomial
logit (179). On simulated data at one thread, with the default scale_vars,
the optimizer's time fell from 14.24 to 7.30 s for a mixed logit (two fixed
and two normal random coefficients, 100,000 choice situations of two to four
alternatives, 50 stored draws), from 8.36 to 4.26 s for a nested logit (two
covariates, 500,000 situations of six alternatives in three nests) and from
4.29 to 3.15 s for a multinomial logit (two covariates, 1,000,000 situations
of two to ten alternatives and an outside option); medians of three
alternated runs. Estimates, standard errors, iteration counts and
convergence codes are unchanged, bit for bit on one thread. nloptr, the
default optimizer, takes the objective and the gradient from one function
and is unchanged.run_mnlogit(), run_mxlogit() and
run_nestlogit(), scaled or not (run_nestlogit()'s default
scale_vars = "none" included), now come from the information matrix
equilibrated before it is inverted, as those of scaled fits already did
(above): scaled on both sides by the powers of two nearest (on a log scale)
the inverse square roots of its diagonal's magnitudes, which adds no
rounding, inverted, and scaled back. This holds at fit time and post hoc
(vcov(fit, type = ), wesml_vcov()), for the inverse of the Hessian or
BHHH information and for the sandwich. The estimator is unchanged. In
natural units an information matrix's curvatures can span many orders of
magnitude (columns in different units, the constants of rarely and often
chosen alternatives); inverted as it was, such a matrix could fail as
computationally singular, leaving the standard errors NA, and a sandwich
built on it could lose digits.NA, with a message, or the reverse.solve(),
with its checks), so an inverse holds two parameters-by-parameters
matrices besides the information where solve() held three, and four or
more with the equilibration. The sandwich is solved twice against the
bread's factors without forming its inverse, and holds two besides the
bread and the meat where it held four, and five with the equilibration.
At 17,000 parameters each such matrix takes 2.3 GB: an inverse holds 4.6
GB besides the information, against 6.9 GB unscaled and 11.6 GB scaled
before, and a sandwich 4.6 GB, against 9.2 and 11.6 GB. On a census-style
mixed logit with 2,620 parameters (9.8 million rows, prepared data), a fit
stopped after one likelihood evaluation peaked at 0.24 GB of heap, its
BHHH variance included, scaled or not, against 0.40 GB unscaled and 0.55
GB scaled before: more than the matrices above account for, since R's own
temporaries of the old inversion stayed in the heap until a garbage
collection. The inverse took the time of solve(), and the
sandwich, at 2,620 parameters, 29% and 38% less at one and four threads
(medians of three alternated runs).run_mnlogit(), run_mxlogit() and run_nestlogit() compute the variance
of their estimates (se_method) once the optimizer has stopped. An error in
that step used to stop the call, and the estimates were lost with it,
however long the optimizer had run: the mixed logit's memory check, say,
which stops a Hessian, BHHH matrix or clustered meat that would need more
than the machine's physical memory (or options(choicer.max_memory =)) at
the threads set, or an allocation that fails in R. The fit now returns its
estimates without a variance (vcov is NULL, the standard errors are
NA, and variance_error holds the error's message) and warns (class
choicer_vcov_failed) with the error and how to compute the variance later:
vcov(fit, type = ), which needs the data the fit keeps with
keep_data = TRUE (the default); the memory check's own message says how
many threads fit, or how to lift its limit. Nothing the size of the variance
is allocated after the failure, so a fit whose variance does not fit in
memory at all keeps its estimates too, as long as the failure is an error
(the memory check, or an allocation that fails with an error). summary(),
print() and the post-estimation methods work as on any fit, with NA
wherever a standard error enters, and do not run the failed step again;
summary() says that the standard errors were not computed and how to
compute them. An allocation that fails while the information is inverted,
which was reported as a singular information (and recomputed at every
summary()), now takes the same route at fit time; after the fit it is an
error, from vcov(fit, type = ) or from summary() and the other methods
when they compute a missing variance. Errors in the optimizer and before it
stop the call as before, and a singular information matrix still gives NA
standard errors with a message. Fits whose variance can be computed are
unchanged, bit for bit. An allocation that fails inside the multinomial and
nested logits' Hessian and BHHH loops, which allocate per thread inside
their parallel region, still terminates R.input_data
that carries choice-based-sampling provenance with uniform weights) now
stops the call before the fit and whose warnings (non-uniform weights
without se_method = "sandwich" or "cluster"; that provenance with
uniform weights) are printed as they occur, before it (in a loop of fits,
once per fit rather than collected at the end); then the fit-time cluster
labels (one per choice situation, none missing, and, for a panel mixed
logit, constant within each decision maker). The messages are the same, and
a call that fails later now gives the guard's warnings as well. Only labels
edited in a prepared input_data can fail the label checks: the
preparations drop the choice situations with a missing value in any column
they use, and check that the labels are constant within each choice
situation and, with person_col, within each decision maker.run_mxlogit() gains polish and polish_control, its last arguments.
With polish = "newton" the optimizer is followed by Newton steps on the
analytic Hessian of the simulated log-likelihood, and with "bhhh" on the
BHHH (outer product of the gradients) matrix, both with the fit's own draws
and weights, from wherever the optimizer stopped. The estimator does not
change: the polish finishes the climb to the maximum the optimizer
approached, which L-BFGS reaches slowly where curvatures differ by orders of
magnitude (thousands of constants, say), and certifies it. It converges when
the Newton decrement nu = g'C^-1 g / 2 (g the gradient, C the curvature)
falls below polish_control$tol_decrement, 1e-8 by default. On a positive
definite Hessian every linear combination of the estimates then lies within
sqrt(2 nu) Hessian standard errors (1.4e-4 at the default) of the maximum of
the log-likelihood's quadratic model, whatever the sample size. BHHH is
positive semi-definite everywhere and cheaper per step, but converges
linearly, and its decrement on its own shows only a stationary point. The
default, "none", leaves every fit as it was, bit for bit.convergence to the optimizer's success code (1 for
nloptr, 0 for optim and a custom optimizer) and leads message with its
verdict. fit$optimizer$polish records the phase (steps, decrement,
damping, evaluations, the parameters on a bound or held, the log-likelihood
gained, its run time) and the optimizer's own code, message, iterations and
run time. A polish that stops short keeps the optimizer's code and says why;
one that stops with an error (the Hessian's memory check, say) keeps its
last accepted point and warns (class choicer_polish_failed). The standard
errors are computed at the polished estimates, by se_method, as before.lower
and upper (projected Newton, Bertsekas 1982), and they pass a line search
that, near the maximum of a likelihood summed over millions of units, also
accepts a step whose gain is below the objective's rounding by reading the
gradient (Hager and Zhang 2005). Each step is solved in one working copy of
the information: at 17,000 parameters the step's peak heap besides the
information was 2.3 GB, against 4.6 GB for a scaled copy factored by
chol(), in about the same time.?run_mxlogit, and §5.7 of the mixed logit's math
article) states the limits. The constants of inside alternatives that have
no finite estimate (next section) are held at the optimizer's values, and
the certificate concerns the other parameters given them; where the outside
option or the base alternative is the one without, the polish does not
report convergence (unless, the reference never chosen, a chosen
alternative's constant held or bounded above stops the others). A random
coefficient whose variance's maximum is at zero drives its log-Cholesky
diagonal toward -Inf: the Newton polish may certify a point whose diagonal,
and its standard error, are set by tol_decrement, and the BHHH polish,
whose curvature there vanishes faster, holds the diagonal and stops
unconverged. "bhhh" needs at least as many likelihood units as the
parameters it moves.run_mnlogit(), run_mxlogit() and run_nestlogit() now warn (class
choicer_unbounded_constants), before the optimizer runs, when the counts
in alt_mapping show that an alternative's constant has no finite
maximum-likelihood estimate: the alternative is never chosen in the data, or
is chosen in every choice situation that offers it. The score of its
constant then has one sign at every value of the parameters, so the
likelihood keeps rising as the constant falls (or grows): the optimizer
stops it wherever its tolerance does, and its standard error means nothing,
while the other parameters approach the maximum without that alternative (or
without the situations that offer it). Where the outside option or the base
alternative is never chosen, or chosen wherever it is offered, the other
constants diverge together. Other patterns of choices can leave constants
without a finite estimate too (a group of alternatives chosen only in
situations that offer no alternative outside the group, say); the counts do
not reveal them, and the warning does not detect them. The estimates are
unchanged.include_outside_option = TRUE adds another. Naming
another base only moves the divergence of a never-chosen base to its own
constant. In the mixed logit a finite bound in the direction a constant
diverges holds it, and that alternative is not named. In the nested logit
the never-chosen case needs the nest's dissimilarity parameter at most 1.
The warning is printed as it occurs, before a fit that can take hours.The following were found while implementing the panel likelihood above and the kernel and data-preparation work for population-scale data that followed, and independently verified. Each affected released versions unless it says otherwise.
mxl_hessian_parallel() silently dropped any choice situation whose
simulated choice probability, summed over draws, was <= 1e-12, while the
objective and gradient kept it. A situation that was merely very poorly fit
by the current parameter vector — not one with genuinely zero simulated
probability — could therefore be silently excluded from the analytical
Hessian, which could then be substantially wrong, and so could the
inverse-Hessian and sandwich standard errors built from it. Such situations
are now kept; only units whose log-likelihood is not finite, which now
takes utilities that overflow (see the log-probability item below), are
skipped.
When a choice situation's simulated probability was subnormal (utility gaps
of several hundred), the gradient came back identically zero alongside a
finite objective value — a false stationary point that could stop the
optimizer early — and the score-based variances (BHHH, robust, cluster)
were non-finite. The posterior draw weights are now normalized in log space
(matching the panel likelihood's omega_ns = exp(lambda_ns - LSE_s lambda_n.)), so these quantities are finite whenever the simulated
probability is positive (and, with the next item, accurate).
The mixed logit kernels took each choice situation's log probability as
log(P). A choice probability that underflowed to zero, at a utility gap
beyond about 745, made its decision maker's log-likelihood -Inf and with
it the whole objective the optimizer's 1e10 sentinel with a zero
gradient; subnormal probabilities, at gaps beyond about 708, lost
precision. The chance of hitting this grows with the number of choice
situations: on a synthetic claims-style panel of 8.7 million stacked rows
with a single choice made against a covariate value of 2000, the fit
returned the sentinel at 342 of 540 evaluations and stopped where that
choice's probability underflows (a coefficient of -0.37, truth -0.5). The
log probability is now the chosen utility minus the log-sum-exp of the
choice set, as in the multinomial logit kernels, so the objective,
gradient, Hessian, scores and conditional tastes are exact whenever the
utilities are finite, however improbable the observed choices; the same
fit now converges (-0.48). The sentinel, the Hessian's skip of a unit and
NA conditional tastes now occur only when the utilities themselves
overflow (for instance an exploding Cholesky factor during a line search),
and because finite objectives are no longer bounded below the sentinel,
run_mxlogit() reports it as ten times the largest objective seen along
the optimizer's path. An explicit kernel overflow flag distinguishes the
sentinel from a valid objective equal to 1e10. Fits away from such
regions agree with earlier versions to optimizer tolerance rather than
bit for bit, since the optimizer's path shifts slightly.
mxl_hessian_parallel() now centers the draw scores in the Louis
identity, sum_s omega_s (g_s - g_bar)(g_s - g_bar)' in place of
sum_s omega_s g_s g_s' - g_bar g_bar'. The two agree in exact arithmetic,
but the former draw weights summed to one only to the precision of their
log-scale arguments, and the uncentered form multiplied that error by
|g|^2: for a choice far below its competitors, an entry of about 1 came
out 0.15 off.
Posterior weights now divide max-shifted exponentials by their sum,
including in the batched score, and Hessian centering first removes a
highest-weight draw's score. This also prevents spurious curvature when a
finite utility gap is so large that subtracting the log-sum-exp loses
normalization accuracy.
The weights described here were this development version's log-space
weights; released versions normalized linear weights, which sum to one to
rounding precision, so there the uncentered form lost accuracy only at far
larger scores.
run_mxlogit()'s advanced workflow (input_data + eta_draws) recorded
the S argument (default 100) instead of the number of draws in
eta_draws, so post-hoc recomputation (vcov(fit, type = ), lazily
computed standard errors, prediction) could regenerate a different draw
set from the one used in estimation. draws_info$S now records
dim(eta_draws)[2]. Post-hoc methods reproduce the estimation draws only
when eta_draws was built with get_halton_normals(S, U, K_w).
The data preparations (prepare_mnl_data(), prepare_mxl_data(),
prepare_nl_data(), prepare_mnp_data(), prepare_hmnl_data(),
prepare_hmnp_data(), and the run_*() functions that call them) could
corrupt or misread an input column whose name matched one of their own.
They added working columns with fixed names (alt_int, idx_in_group,
HAS_NA, and in the hierarchical preparations HB_PERSON, task_idx,
TASK_HAS_NA, HAS_BAD, TASK_HAS_BAD) to a copy of the data, read
per-situation results back under names (N, V1, chosen, pos) that
an id column of the same name shadowed, and looked up variables of their
own (levels, J, outside_opt_label, ids_to_drop, id_col,
choice_col, and in the hierarchical preparations person_col and
task_by) where a column of the same name took their place. predict(),
logsum() and consumer_surplus() with a data.frame newdata,
wesml_weights(), sample_by_choice(), and the hierarchical fits and
their post-estimation methods did the same: they looked up variables of
their own (am, pos, spec, d, choice_col, id_col, keep_ids)
inside the data or the fit's alt_mapping, and read counts and strata
back under names (V1, .strat) that an id column of the same name
shadowed. These clashes changed the results without an error:
prepare_hmnl_data() / prepare_hmnp_data(), a covariate or
control-function residual named alt_int, idx_in_group or
task_idx, or one named HB_PERSON when the respondent ids were
numeric, was replaced by the alternative codes, the within-situation
positions, the task indices or the respondent ids (values -1.2, -0.69,
2.29 entered the design matrix as 1, 2, 3), as was an alternative-level
covariate named idx_in_group when every choice situation offered the
same alternatives; prepare_mnp_data() differenced the codes of a
covariate named alt_int;task_by, other than the
task column of a fit without person_col, regrouped the rows by its own
values: a covariate or control-function residual whose values differ
from row to row made every row a choice situation of its own, without
an error whenever the outside option was modelled, as it is by default
(repeated values usually stopped the preparation);prepare_mnl_data(), prepare_mxl_data() and prepare_nl_data(),
a choice column also listed in covariate_cols left the working table
before the choices were counted: the preparation stopped with an
unrelated error or, when an object of the column's name was visible to
the call (a vector choice in the workspace, say), counted that object
instead, so that alt_mapping credited every alternative with every
choice and, unless the rows were already in prepared order, the chosen
alternatives were wrong;alt_int or idx_in_group recorded the first
alternative of every choice situation as chosen (idx_in_group was
harmless in prepare_mnp_data()), and in the hierarchical preparations
one named task_idx recorded the outside good in every task but the
first;alt_int or idx_in_group with the outside option,
and in the hierarchical preparations a task column named alt_int or
task_idx, regrouped the rows into spurious choice situations;outside_opt_label kept the outside-option
rows as an extra inside alternative;idx_in_group with unequal choice sets
added unidentified constants, leaving every standard error NA;predict(), logsum() and consumer_surplus() with a data.frame
newdata, a covariate named pos whose values were valid positions in
the fit's alt_mapping (without an outside option, any permutation of 1
to J within each choice situation, such as a rank) gave each row the
alternative at that position, so the alternative-specific constants and
the shares went to the wrong alternatives;keep_ids made sample_by_choice() return the wrong
choice situations (the whole population when it was the id column).These changed labels and summaries: an alternative column named
alt_int, or with equal choice sets idx_in_group, had its labels
replaced by the codes; one named TAKE_RATE or MKT_SHARE had its
labels, and so the ASC names, overwritten by the take rates or shares;
one named N_OBS or N_CHOICES with numeric labels misstated the take
rates or shares, and made gof(null = "market_shares") stop (N_OBS)
or return a wrong null log-likelihood (N_CHOICES); and one named J
mislabelled the constants of run_mnlogit(), run_mxlogit(),
run_nestlogit() and prepare_mnp_data() (with three alternatives
labelled 1 to 3 or stored as a factor, the third alternative's constant
took the base alternative's label). Whatever matched alternatives by
label inherited the error: prepare_nl_data() assigns nests by label, so
it stopped or, with text labels such as "1" to "10", assigned wrong
nests; predict(), logsum() and consumer_surplus() with newdata
stopped or attached another alternative's constant (labels 0 to 2 had
become 1 to 3); and predict() on the hierarchical fits, which takes an
unknown label for a new alternative, predicted every alternative from
posterior-predictive constants without an error.
An id column named N gave the ids as the choice-set sizes M, in the
preparations and in predict(), logsum() and consumer_surplus() with
newdata. The kernels' checks usually stopped the fit or prediction, but
a prediction whose ids summed to the number of rows (ids 1 to 11 with 6
alternatives each, for instance) ran on wrong choice sets without an
error, and more rarely so could a fit. A covariate named levels in
prepare_mnl_data(), prepare_mxl_data() and prepare_nl_data() left
every alternative code NA (or, when its values repeated, stopped the
preparation), and with the outside option an id column named pos gave
the ids as the chosen positions; fitting then failed.
Other clashes stopped with unrelated errors: an id column named V1 (the
name data.table::fread(header = FALSE) and as.data.frame() of an
unnamed matrix give the first column) in prepare_mnl_data(),
prepare_mxl_data() and prepare_nl_data() whenever a weight or cluster
column was given, and always in the hierarchical preparations; a
covariate named alt_int or idx_in_group in prepare_mnl_data(),
prepare_mxl_data() and prepare_nl_data(); a column named
choice_col that a model used (other than the nested logit's nest
column), as did one named id_col outside the hierarchical preparations
and, in those, one named person_col when person_col was given or
id_col when it was not; an id or alternative
column whose name contains a comma, which data.table's grouping splits
into several names; and an id, weight or cluster column named like a
working column, among others. So did, outside the preparations, a
covariate named am in a data.frame newdata, most covariates there
named pos, and with an outside option one named spec; an id column
named V1 or .strat, or a choice column named choice_col, in
wesml_weights() and sample_by_choice(); and an alternative column
named am, which stopped run_hmnlogit() and run_hmnprobit(), or
d, which stopped their predict() without newdata, elasticities(),
diversion_ratios() and ppc_shares(), and for the HMNL logsum()
without newdata and consumer_surplus(). Fits whose data used a name
in the first group, and fits with an id column named N that ran, should
be re-estimated, and the predictions and samples in that group
recomputed; for the second group, re-estimate nested logit fits and
recompute predictions made with newdata.
Covariates are now read from the data rather than from the working table;
the remaining working columns, per-situation results, counts and strata
carry the reserved prefix .choicer_ or are read by position; and the
preparations, post-estimation and the WESML helpers compute their own
values, and read the columns they are given by name, before indexing the
data or alt_mapping, so the names listed above are no longer looked up
among the columns. Inputs without a clash give the same results as
before. A choice column stored as text or as a factor, which the
preparations, wesml_weights() and sample_by_choice() never accepted,
and in prepare_mnl_data(), prepare_mxl_data() and prepare_nl_data()
a covariate that is also the id, alternative, choice, weight, cluster or
decision-maker column (see above for the choice column), now stop them
with an error that says why. Two inputs are now errors: a column whose
name starts with .choicer_, used as an id, alternative, choice,
covariate, weight, cluster, nest or decision-maker column; and an
alternative column named alt_int, N_OBS, N_CHOICES, TAKE_RATE or
MKT_SHARE, the fixed columns of the returned alt_mapping, which must
now be renamed. Released versions handled .choicer_ columns correctly
(they used no such names) and, of the alternative-column names, only an
alt_int column already coded 1 to J, without the outside option.
A covariate named like a parameter the fit generates gave two parameters
one name. prepare_mnp_data() (and so run_mnprobit()) names the
constant of each non-base alternative ASC_<label> and builds param_map
by column name: since 0.2.0, a covariate such as ASC_b next to
alternative b's constant left that constant in neither the beta nor
the asc block, so recovery_table() could not align the ASCs, and
summary() stopped with a duplicate row-name error; when the covariate
was itself alternative b's dummy, one of the two columns was dropped as
collinear without the usual message. run_mnlogit(), run_mxlogit() and
run_nestlogit() generate ASC_<label>, Lambda_<k> (nested logit),
Mu_<variable> and L_<i><j> (mixed logit), and the mixed logit's
summary() prints Sigma_<i><j> for L_<i><j> and exp(Mu_<variable>)
for the mean of a log-normal coefficient; a covariate could take any of
these names. Estimation addresses their parameters by position, so
without named bounds (below) the estimates, standard errors and
predictions were right. But summary() stopped with a duplicate
row-name error, unless it relabelled the generated parameter
(L_<i><j>, or a log-normal Mu_<variable>),
and indexing coef() or vcov() by name, or wtp(attr_vars = ), found
the covariate: wtp(attr_vars = "ASC_b") returned the covariate's WTP,
and the constant's could not be reached by name. A named bound in
run_mxlogit() goes to the first parameter of that name, and the fixed
coefficients come first, so it bounded the covariate's coefficient and
left the generated parameter free. With a covariate named ASC_b,
lower = c(ASC_b = 0) held that covariate's coefficient at 0 whenever its
unconstrained estimate was negative, and so moved the other estimates,
while alternative b's constant stayed unbounded; a covariate named
L_11 would likewise have taken the documented lower = c(L_11 = -5)
meant for the first Cholesky diagonal. Such names are now an error that
lists them: prepare_mnp_data() stops before building X (rename the
covariates, or set use_asc = FALSE for ASC dummies built by hand), and
the three other functions stop before optimizing, run_mxlogit() before
building its draws. The same check stops an input_data whose X
repeats a column name, and repeated run_nestlogit(param_names = ), which
used to fit with the same ambiguities. Preparations and fits whose
parameter names and summary() labels are all distinct are unchanged,
summary() included.
Columns of class integer64 (bit64), which data.table::fread() returns
for integers beyond 2^31 - 1 and several database drivers return for
BIGINT, passed the check that covariates are numeric but entered the
design matrices as their raw 64-bit patterns read as doubles, in every
model of released versions. Positive values became tiny numbers (1 read as
4.9e-324, 2^40 as 5.4e-312), so fits, standard errors and predictions were
silently wrong. Negative values above -2^52 became NaN, which stopped
the preparation with an unrelated error (NA/NaN/Inf in foreign function call); with bit64 loaded, prepare_mnp_data() differenced these columns
as integers and so usually stopped the same way. When bit64 was not
loaded, as for data read from an .rds file, the check for missing
values misread these columns too: it took negative values for missing
ones, dropping their choice situations with the warning about missing
values, and let a missing value into the design as 0. The design
builders now read integer64 columns as their values, as
as.double() converts them: prepare_mnl_data(), prepare_mxl_data(),
prepare_nl_data(), prepare_mnp_data(), prepare_hmnl_data() and
prepare_hmnp_data() (including alt_covariate_cols and
cf_residual_col), and the newdata of predict(), logsum() and
consumer_surplus(). The conversion is exact below 2^53 in magnitude;
from 2^53 up, values are rounded to the nearest double, with a warning
naming the column. An integer64 column among those a model uses now loads
bit64 first, so its missing values are found and handled as in other
columns; if bit64 is not installed, such a column is an error. Weights
(weights_col, a weights vector) and integer64 X and W matrices in
the list form of newdata had the same defect and are now converted the
same way, as are prediction weights. bit64 is now a suggested package.
sample_by_choice() misread integer64 ids, as data.table::fread()
returns ids beyond 2^31 - 1, ever since it was added in 0.2.0: it matched
the sampled ids back to the rows after unlist() had dropped their
class. With bit64 loaded, the sample came back empty without an error;
without it, each id drawn from among ids of large magnitude also selected
its neighbours (about six choice situations for each one requested, for
ids near 10^18). With
the outside option, wesml_weights() returned its ids as raw bits
(doubles such as 7.8e-242) whenever a choice situation chose the outside
good and bit64 was not loaded, and with attach = TRUE stopped with an
unrelated error. Both functions now load bit64 when the id, alternative
or choice column is integer64, as the preparations do, and keep the
ids' class, so they sample and weight such ids by value; if bit64 is not
installed, such a column is an error. Samples of other ids are drawn as
before.
nl_blp_contraction() takes lambda and nest_idx as given, not from a
fitted model, and never checked them against each other: a nest code
below 1 or missing, or a lambda with other than one entry per nest,
failed an Armadillo index or size check inside the contraction's parallel
loop, which aborted R in builds with OpenMP and was an Armadillo error in
builds without it. Both are now errors before the contraction starts. As
in the other nested logit kernels, this includes a nest code below 1 on
an alternative that no choice situation offers, which used to be
ignored.
run_mnprobit(), run_hmnlogit() and run_hmnprobit() miscounted the
draws to keep when R - burn + thin - 1 exceeded .Machine$integer.max,
as for a thin near it: the count overflowed an integer. Where it came
out 0 this aborted R in builds with OpenMP; otherwise the run stopped
with an unrelated Armadillo error (with keep_beta_i = "draws", the
memory estimate failed first, with another unrelated error). The count,
ceiling((R - burn) / thin), is now formed without overflow.
Mixed logit fits with stored draws (draws = "store", the default) read
some units' draws from the wrong place, without an error, once the draw
cube (K_w * S * N values) was large enough that a unit's draws began
past its first 2^32 values. The cube was built for any S * N up to
2^31 - 1 points (past that, building it failed). With three random
coefficients, the default 100 draws and 15 million choice situations, for
instance, it holds 36 GB, and the last 683,442 situations read other
units' draws in place of their own. Those draws entered the estimates and
their standard errors, and every result that rebuilds the cube: standard
errors recomputed with vcov(type = ) or wesml_vcov(), predictions,
elasticities, diversion ratios, BLP inversions, log-sums, consumer surplus
and fit statistics. Such results should be recomputed. The compiled code
used RcppArmadillo's default 32-bit word, under which the kernels'
in-place view of a matrix or cube of more than 2^32 - 1 values wrapped its
element count and element offsets: a unit's draws, or a choice situation's
covariate values, that began past the first 2^32 values were read from the
wrong place. A hand-built eta_draws that large was misread the same way.
A design matrix that large (for instance 10^9 rows and 5 covariates)
reached the kernels only through a hand-built input_data, the newdata
of a prediction or welfare function, or a direct call to an exported
contraction function: the preparations stopped it earlier, at qr()'s
limit of 2^31 - 1 values. The compiled code now uses Armadillo's 64-bit
word (ARMA_64BIT_WORD), so the kernels read matrices and cubes of any
size R can hold, and outputs of more than 2^32 - 1 values, such as the
output of conditional_tastes(), are allocated where Armadillo stopped
with "requested size is too large". Below 2^32 values, results are
unchanged, bit for bit.
summary() of a mixed logit fit without a variance, one whose information
matrix is singular, printed the Cholesky factor's parameters (the diagonal
on its log scale) under the labels of the covariance matrix,
Sigma_<i><j>, and the means of log-normal coefficients untransformed
under exp(Mu_<variable>): it transformed the estimates only when it had a
variance for their standard errors. It now reports the covariance matrix and
exp(Mu) whether or not the variance is available, with NA standard
errors where it is not.
run_mxlogit(use_asc = FALSE), the mixed logit without alternative-specific
constants, counted the J - 1 constants of the model with them whatever
use_asc said. Since 0.2.0 the first evaluation of the likelihood, at start
values that many entries too long, stopped with "Theta vector too long: ...
Did you mean use_asc = TRUE?". 0.1.0, whose kernels did not check the
length, returned the constants as coefficients held at their start of 0, on
which the likelihood does not depend: its other estimates were those of the
model without constants, but n_params, AIC and BIC counted the constants,
and every standard error was NA (the information is zero in their rows and
columns, so singular). It now fits the model without constants: its
parameters are the fixed coefficients, the random coefficients' means (with
rc_mean = TRUE) and their Cholesky factor, param_map has no asc block,
as in run_mnlogit() and run_nestlogit(), and the default start,
theta_init and lower/upper have one entry per parameter, with or
without the outside option. Prediction, elasticities, diversion ratios, BLP
inversions, log-sums, consumer surplus, WTP, conditional tastes and every
vcov(type = ) work on such fits. Fits with constants (use_asc = TRUE,
the default) are unchanged, bit for bit.
blp() on a model fitted without constants (use_asc = FALSE) and with
the outside option (include_outside_option = TRUE) started its
contraction from one zero too many, counting the outside option, which
has no constant, among the alternatives. The multinomial and nested logit
methods stopped with an error that blamed target_shares ("must have the
same length as the total number of alternatives") ever since blp()
gained them (the multinomial logit in 0.1.0, the nested logit in 0.2.0);
the mixed logit's had the same start, which no fit of a released version
reached. The contraction now starts from one zero per inside alternative;
other fits are unchanged.
predict(type = "shares") with newdata that does not offer the
alternatives that come last in object$alt_mapping, as when an
alternative exits in a counterfactual, returned one share fewer per such
alternative for multinomial and nested logit fits, with or without
constants, since newdata was added in 0.2.0; mixed logit fits without
constants, possible since the use_asc = FALSE correction above, did the
same. The shares were right, but the result stopped at the last
alternative that newdata offers, so it no longer lined up with the
baseline's shares; a mixed logit fit with constants, and any model when an
earlier alternative is missing, gave each alternative that newdata does
not offer a 0. Multinomial, nested and mixed logit fits now return one
share per alternative of object$alt_mapping, in its order (the outside
option's first, when it is modelled), 0 for an alternative that newdata
does not offer; other predictions are unchanged, bit for bit. The
documentation called the result a named numeric vector; it is, as it
always was, a one-column matrix without names (the labels are in
object$alt_mapping[[2]]).
An outside_opt_label that names none of the alternatives, without the
outside option (include_outside_option = FALSE), became a phantom base:
the preparations code that label first, so it took code 1 with no rows.
Fits with constants then stopped with an error about the constants' block
(with three alternatives, "Theta's delta (ASC) block implies 3
alternatives but alt_idx references alternative 4"; an Armadillo index
error in 0.1.0), nested logit fits without them with an error about the
nests (in 0.1.0 R aborted), and multinomial logit fits without constants,
like mixed logit fits without them since the use_asc = FALSE correction
above, carried an alternative that is never offered: a predicted share of
0, and elasticity and diversion-ratio matrices with an extra row and
column and no labels. An NA label silently made the first alternative
the base. A label of another class than the alternative column (a
character label with a column of dates, a plain number with an
integer64 column, a factor label with a character column) could leave
alternatives without a code, which the preparations returned: with dates
and a character label, checked on Apple silicon, fits since 0.2.0 stopped
with an error about the alternative codes, and 0.1.0 fitted a model
without constants on them. prepare_mnl_data(), prepare_nl_data() and
prepare_mxl_data(), and the run_*() functions that call them, now stop
in these cases with an error that says why. Labels that coded the
alternatives before still do: the check reads the codes the preparations
build.
Patch release addressing a compilation warning reported by CRAN's GCC check flavors. No user-visible behavior, API, or numerical results change.
#pragma omp master directives in the Gibbs samplers and
thread-introspection helpers with a CHOICER_OMP_MASKED macro that emits
masked on OpenMP 5.1 and later, master on earlier versions, and nothing
when the package is built without OpenMP. OpenMP 5.1 deprecated master in
favor of masked, and GCC 16 — now the compiler on the
r-devel-linux-x86_64-fedora-gcc and r-devel-linux-x86_64-debian-gcc
check flavors — warns on the old spelling, which R CMD check reports as a
significant install warning, escalating the check result to WARNING.
masked with no filter clause is semantically identical to master: the
block runs on the primary thread with no implied barrier, so the samplers'
fixed-order reductions, barrier structure, and
no-R-API-off-the-primary-thread contract are unchanged, as are posterior
draws for a given seed and thread count.get_halton_normals(),
set_num_threads(), and thread_info() remain public. With no downstream
CRAN dependencies, v0.2.0 is the least disruptive point to narrow this
surface before applications depend on unstable kernel signatures.vcov() on a fitted MNL, MXL, or NL model gains type = and cluster =
arguments for post-hoc variance recomputation without refitting (requires
keep_data = TRUE): "hessian" (inverse analytical Hessian), "bhhh"
(OPG), "robust" (Huber-White / WESML sandwich), and "cluster"
(cluster-robust sandwich over within-cluster sums of weighted scores — use
it when the same decision maker contributes several choice situations).
With no arguments, vcov() returns the as-fitted variance, unchanged.cluster_col= argument on run_mnlogit() / run_mxlogit() /
run_nestlogit() (and the corresponding prepare_*_data() functions):
supplies per-situation cluster labels at fit time and selects the new
se_method = "cluster".se_method = "sandwich" /
wesml_vcov() results are unchanged (the robust meat is the w^2 special
case of the shared path).run_hmnlogit() (person_col=).blp() for an MXL fit using draws = "generate" now carries the fitted
draw count, seed, and digit-permutation setting through every contraction
evaluation. Previously that path could evaluate counterfactual shares with
fallback draw metadata rather than the simulation design used for the fit.person_col, id_col) before subtraction;
it now rejects added, dropped, or substituted tasks instead of pairing them
silently by sorted position.run_mnprobit(), run_hmnlogit(), run_hmnprobit()) are
bitwise reproducible regardless of the OpenMP thread count. Independent
validation disconfirmed this: draws are reproducible given the seed and a
fixed thread count, and across different thread counts are invariant only
up to floating-point reduction-order round-off (~1e-15), not bitwise. All
man pages, vignettes, and NEWS entries now state the correct guarantee.run_mxlogit() as a
cross-sectional simulated likelihood, and no longer directs users to
WTP-space estimation or bounded/censored mixing distributions as if choicer
implemented them. The README also labels mode_choice shares and welfare as
sample-design quantities unless external population shares and WESML weights
are supplied.run_mnprobit().ess(draws) (rank-normalized bulk and tail effective sample size) and
mcse(draws, kind = c("mean", "median")) (Monte Carlo standard error of
the mean or median), following Vehtari, Gelman, Simpson, Carpenter &
Bürkner (2021, Bayesian Analysis). rhat() gains a rank = TRUE option
for the same paper's rank-normalized, folded R-hat; the default
rank = FALSE reproduces the original split R-hat exactly.traceplot() generic, with a choicer_hb method, for visual
inspection of the b/theta/sigma_d2 chains (and user-selected delta
columns).summary() on a choicer_hmnl/choicer_hmnp fit now prints one
consolidated multi-chain diagnostics table (R-hat, ESS bulk, ESS tail, and
MCSE per parameter block, plus a single worst-case summary row spanning
all J delta_j alternative effects) in place of the previous bare
split-R-hat printout. HMNP now always prints an explicit
"conjugate — no acceptance step" line where HMNL prints its beta/delta
acceptance rates. The fit-time convergence warning now checks rank-R-hat
and ESS bulk across every tracked parameter, including all J delta_j
columns (previously only b/theta/sigma_d2 were checked).run_hmnlogit() / run_hmnprobit() fits now retain all requested
chains' hierarchical draws in a new object$chains field (a list, one
element per chain: b, w_vech, delta, theta, sigma_d2, plus
loglik for HMNL / sigma2 for HMNP), which the multi-chain diagnostics
above consume. object$draws (chain 1) and object$beta_i
(chain-1-only respondent taste summaries/draws) are unchanged.chains= are sampled sequentially, and all of their
hierarchical draws are retained (object$chains) for the multi-chain
diagnostics above. Parallel multi-chain execution is planned for a future
release.keep_beta_i = "draws" memory guard on both models is now based on a
measured ~1.9x Rcpp::wrap()/R-list overhead factor (the previous formula
under-estimated true memory use by roughly 2x) and accounts for all
requested chains; it still fails fast before the chain runs rather than
after an out-of-memory blow-up.run_hmnlogit() (class choicer_hmnl) and run_hmnprobit() (class
choicer_hmnp): hierarchical Bayes discrete choice with two decoupled
random-effect levels — respondent-level structural tastes
beta_i ~ N(b, W) (W is K x K, no ASC columns) and a global BLP-style
alternative effect delta_j = z_j' theta + xi_j, xi_j ~ N(0, sigma_d^2)
(one scalar variance regardless of J; partial pooling toward the
characteristics-based mean; posterior-predictive delta for alternatives
outside the estimation sample). Both models carry a first-class implicit
outside option (systematic utility 0 plus its own shock) that anchors the
location of delta — no base alternative or sum-to-zero constraint.
Panel and cross-sectional (person_col = NULL, T_i = 1) modes share
one code path.delta sweep
(the conditionals are coupled through the softmax denominators) at O(1)
incremental cost per affected task; supports log-normal coordinates via
rc_dist. HMNP: fully conjugate Albert-Chib augmentation in
un-differenced utility space with iid N(0, sigma^2) shocks, a
non-identified sigma^2 chain (parameter expansion), and per-draw scale
normalization of every identified quantity.b ~ N(b_bar, A^-1), W ~ IW(nu, V),
theta ~ N(theta_bar, A_theta^-1), and sigma_d ~ half-Cauchy(0, s_d)
via the Makalic-Schmidt conjugate scale mixture (IG fallback available).prepare_hmnl_data() / prepare_hmnp_data() (two-level
person/task indexing, implicit-outside convention shared with
prepare_mnl_data(), alternative-level design Z, optional
control-function residual column for price endogeneity per Petrin &
Train 2010) and DGPs simulate_hmnl_data() / simulate_hmnp_data().choicer_hb class, all posterior-draw
based: predict() (population/individual, entry counterfactuals via the
posterior-predictive delta; HMNP probabilities by deterministic 20-node
1-D Gauss-Hermite approximation), wtp() (posterior median + quantile
intervals), logsum() / consumer_surplus() (HMNL-only; probit Emax is
roadmapped), elasticities() / diversion_ratios() (common-random-path
perturbation engine including the outside option), recovery_table(),
plus rhat() (split R-hat) and ppc_shares() diagnostics.vignettes/articles/hierarchical_mnl_math.Rmd
and hierarchical_mnp_math.Rmd; recovery demos in
inst/simulations/hmnl_simulation.R / hmnp_simulation.R.run_mnlogit() and run_nestlogit() now compute the robust (Huber–White /
WESML) sandwich variance V = A^{-1} B A^{-1}, at full feature parity with
run_mxlogit(). Both gain se_method = "sandwich" (bread = weighted negated
Hessian, meat = weight-squared OPG) and se_method = "bhhh" (ordinary
outer-product-of-gradients). run_mnlogit()'s se_method therefore now
accepts "hessian" (default), "bhhh", or "sandwich"; run_nestlogit()'s
accepts "hessian" (default), "numeric", "bhhh", or "sandwich".mnl_bhhh_parallel() and nl_bhhh_parallel() accumulate the
weighted outer product of per-individual scores. Their per-individual score is
weight-free, so passing weights = w yields the BHHH information and
weights = w^2 yields the sandwich meat. The NL kernel includes the full
beta/lambda/delta score blocks (singleton-nest lambdas fixed to 1 contribute
no score).run_mnlogit(), run_nestlogit(), prepare_mnl_data() and
prepare_nl_data() gain a weights_col argument: a row-level weight column
is collapsed to one weight per choice situation (validated constant within
id). A choice-based-sampling provenance guard auto-adopts the recorded
WESML weight column and errors rather than silently fitting unweighted under a
WESML label; non-uniform weights under a non-sandwich se_method emit a
warning.se_method = "sandwich" when non-uniform weights are supplied with the default
("hessian") or "bhhh" method; the "bhhh" case gets a sharper message
explaining that BHHH/OPG is not a valid WESML correction (its meat is w^1,
not w^2). Point estimates and standard errors are unchanged.wesml_vcov() now dispatches on choicer_mnl and choicer_nl (in addition
to choicer_mxl), returning the post-hoc sandwich variance from a fit stored
with keep_data = TRUE.summary() for MNL and NL fits now reports the standard-error method and any
WESML weighting in the printed footer.prepare_mnl_data(), prepare_nl_data(), and prepare_mxl_data() (covering
both the weights= and weights_col= paths). Zero, negative, or non-finite
weights previously could slip through and silently invalidate weighted and
WESML sandwich inference (weight w enters the bread, w^2 the meat); they
now error with an actionable message.input_data= passed directly to run_mnlogit(),
run_nestlogit(), or run_mxlogit()) that carry WESML choice_sampling
provenance but resolve to uniform weights now error instead of warning. The
message explains how to proceed: bake the non-uniform WESML weights into
input_data via prepare_*_data(weights=/weights_col=), or strip the
provenance with attr(input_data, "choice_sampling") <- NULL for a deliberate
unweighted fit. Convenience-mode behavior is unchanged.mode_choice data set: the classic Greene & Hensher intercity
travel-mode choice data (210 travellers x 4 modes), in choicer's long layout,
used by the getting-started vignette. ?mode_choice and the vignettes now
document that the sample is choice-based (car under-sampled): slopes and WTP
ratios are unaffected in the ASC-saturated logit, while constants, shares,
and surplus levels inherit the design — see the WESML vignette for the
correction.blp() as the Berry 1994
inversion), what identifies the nested-logit dissimilarity parameters, the
Keane (1992) covariance-identification caveat for the multinomial probit,
and the classical MSL asymptotics (fixed-S bias, S growth conditions) in
the mixed logit math note.run_mxlogit() gains three new arguments: draws, seed, and scramble.
draws = "store" (default) keeps the existing behavior: a full
K_w × S × N Halton cube is pre-materialized and stored in memory.draws = "generate" activates the new mode: each individual's S draws are
computed on the fly in C++ from a compact seed, eliminating the O(N)
eta_draws cube. This is the recommended choice when N is large or memory
is constrained.scramble = "permuted" (default when draws = "generate") applies one
seeded permutation per dimension and base-digit position, shared across
sequence indices. This is a deterministic position-wise digit permutation,
not Owen's nested-uniform scramble, and no standard randomized-QMC
unbiasedness or replicate-error guarantee is claimed. The historical value
"owen" remains as a deprecated alias for compatibility.
scramble = "none" reproduces the randtoolbox sequence exactly.seed sets the integer master seed for the on-the-fly generator; if
NULL (default), a seed is drawn from R's RNG so set.seed() governs
reproducibility. Ignored when draws = "store".draws_info field on choicer_mxl objects gains three new elements:
mode ("store" or "generate"), seed, and scramble. Existing code
that accesses draws_info$S, draws_info$N, or draws_info$K_w is
unaffected; the new fields are NULL for objects fitted before this release.predict(), elasticities(),
diversion_ratios(), logsum(), consumer_surplus(), vcov()) propagate
the fitted draw mode automatically — no user action required.draws = "store") is bitwise unchanged.run_mnprobit() — Bayesian multinomial probit via Gibbs sampling with data augmentation (Albert & Chib 1993; McCulloch & Rossi 1994). Runs the non-identified chain with conjugate priors and reports identified quantities normalized per draw by sigma_11. New choicer_mnp posterior object with summary() (posterior mean, SD, credible intervals), coef(), vcov(), nobs(); math note in vignettes/articles/bayesian_multinomial_probit_math.Rmdsrc/rng.h, src/bayes_samplers.h): xoshiro256++/splitmix64 RNG with one stream per (iteration, observation) — draws are reproducible given the seed and a fixed OpenMP thread count (invariant across thread counts only up to floating-point reduction-order round-off, ~1e-15) — plus exact truncated-normal, multivariate-normal, Wishart (Bartlett), and inverse-Wishart samplers. The truncated normal picks the cheapest exact method per region (naive normal rejection in high-mass regions, Robert 1995 exponential rejection in the tail, inverse CDF for narrow intervals)_benchmarks/ MNP preset this samples ~2x faster than MNP::mnp() and ~3.5x faster than bayesm::rmnpGibbs() on one thread, and scales with threads on topsimulate_mnp_data() — probit DGP returning a choicer_sim with truth on the identified scale; recovery_table() gains a choicer_mnp method (posterior mean/SD, normal-approximation credible intervals, sigma block); parameter-recovery walkthrough in inst/simulations/mnp_simulation.Rwtp() — willingness-to-pay for MNL, MXL, and NL with analytic delta-method standard errors. For MXL, log-normal random coefficients report the median WTP under the package's shifted log-normal parameterization; random price coefficients are rejectedgof() — goodness of fit: McFadden pseudo R-squared (plain and adjusted) and in-sample hit rate, with "equal_shares" (default) and "market_shares" null models; now also shown in the summary() footerpredict(..., newdata = ) — counterfactual prediction for all three models, from either a long data.frame in the fit-time format or a modified-design list (X, alt_idx, M, ...); works even with keep_data = FALSE. NL fits now store nest_idx top-level to support thislogsum() — expected maximum utility (inclusive value) per choice situation for MNL (closed form), MXL (simulated with a dedicated per-draw kernel, avoiding the Jensen bias of averaging utilities first), and NL (nested inclusive-value formula)consumer_surplus() — expected consumer surplus logsum / (-alpha) (Train 2009, Ch. 3) with a delta-method standard error of the mean CS for MNL; supports newdata for policy ΔCS analysisInitial CRAN release.
run_mnlogit()) — estimation, prediction, elasticities, diversion ratios, BLP contractionrun_mxlogit()) — normal and log-normal random coefficients, correlated random coefficients via Cholesky, Halton draws, elasticities, BLP contractionrun_nestlogit()) — estimation with nest-specific dissimilarity (lambda) parameterschoicer_fit with subclasses choicer_mnl, choicer_mxl, choicer_nlsummary(), coef(), vcov(), logLik(), AIC(), BIC(), nobs(), predict()choicer_data_mnl, choicer_data_mxl, choicer_data_nl from prepare_*_data()elasticities() — methods for MNL and MXLdiversion_ratios() — method for MNLblp() — BLP contraction for MNL and MXLrun_*logit() functions: convenience (pass data + column names) or advanced (pass pre-prepared input_data)optimizer = "nloptr" | "optim" | <function>prepare_nl_data() added for nested logit data preparation