# Bayesian multinomial probit

The multinomial probit (MNP) replaces logit’s Type-I extreme-value error structure with a multivariate-normal covariance for utility differences. That covariance can encode correlated substitution across alternatives. The cost is scale normalization, a covariance parameterization that grows quickly with the number of alternatives, and MCMC diagnostics. choicer estimates the model the Bayesian way: a Gibbs sampler with data augmentation, written in C++ with a reproducible, thread-safe random number generator.

    library(choicer)
    set_num_threads(2)

## Simulate from a probit process

`simulate_mnp_data()` draws choices with correlated normal errors and known parameters.

    sim <- simulate_mnp_data(N = 2000, J = 3, seed = 1)
    sim
    #> <choicer_sim: mnp>
    #>   settings:
    #>     N = 2000
    #>     J = 3
    #>     K_x = 2
    #>     base_alt = 1
    #>   rows in $data: 6000
    #>   true_params: beta, delta, Sigma

## Run the sampler

`run_mnprobit()` returns posterior draws. The settings below keep this vignette quick; for real work use a longer run, repeat the fit with several independent seeds when needed, and inspect trace, autocorrelation, ESS and Monte Carlo error before interpreting posterior summaries. `run_mnprobit()` does not currently expose user-specified starts or a multi-chain wrapper.

    set.seed(3)
    fit <- run_mnprobit(
      data           = sim$data,
      id_col         = "id",
      alt_col        = "alt",
      choice_col     = "choice",
      covariate_cols = c("x1", "x2"),
      mcmc           = list(R = 4000, burn = 1000, thin = 2)
    )
    #> MCMC run time 0h:0m:0.54s
    summary(fit)
    #> Bayesian Multinomial Probit (MNP) model
    #>
    #> Parameter        Mean         SD       2.5%     Median      97.5%
    #> x1           0.784144   0.046153   0.695657   0.783920   0.876807
    #> x2          -0.583644   0.041269  -0.664795  -0.581002  -0.504487
    #> ASC_2        0.472953   0.043467   0.383959   0.473618   0.558991
    #> ASC_3       -0.537857   0.108212  -0.791786  -0.522078  -0.359531
    #>
    #> Covariance of utility differences (Sigma, identified scale):
    #> Parameter        Mean         SD       2.5%     Median      97.5%
    #> Sigma_11     1.000000   0.000000   1.000000   1.000000   1.000000
    #> Sigma_21     0.410913   0.134128   0.119895   0.422711   0.660771
    #> Sigma_22     1.410768   0.277005   0.954727   1.387811   1.997864
    #>
    #> Posterior mean Sigma:
    #>        w_2    w_3
    #> w_2 1.0000 0.4109
    #> w_3 0.4109 1.4108
    #>
    #> Base alternative: 1
    #> Draws kept: 1500 (R = 4000, burn = 1000, thin = 2, seed = 721735354)
    #> N: 2000  | Parameters: 4
    #> Sampling time: 0.54 s
    #> Identification: per-draw normalization by sigma_11 (McCulloch-Rossi 1994).

The `Sigma` entries are the error-covariance parameters. They are identified only up to scale, so choicer reports them on the normalized scale where `Sigma_11 = 1`. The prior is placed on the unrestricted covariance used inside the Gibbs sampler; the reported posterior summaries are computed after normalizing each kept draw.

For empirical work, ask what disciplines the covariance before interpreting it. Keane (1992) shows that when covariates do not vary across alternatives — demographics, or constants alone — the covariance parameters are identified by the normal functional form only, and estimation is fragile. Alternative-specific covariate variation (prices, distances, or attributes that differ across the options within a choice situation), and especially defensible exclusions across utility equations, is what separates covariance from systematic utility in practice; the simulated data above contain two varying attributes. The Bayesian machinery does not repeal this: with weak alternative-specific variation the \(\Sigma\) posterior sits close to its prior, and the substitution pattern it implies is maintained rather than estimated. The MNP is attractive when flexible error correlation is the central object, but with many alternatives its parameter count rises quickly.

## Compare posterior summaries with the truth

    recovery_table(fit, sim$true_params)
    #> <choicer_recovery> model=choicer_mnp level=0.95
    #>    parameter  group  true estimate     se    bias rel_bias_pct z_vs_true
    #>       <char> <char> <num>    <num>  <num>   <num>        <num>     <num>
    #> 1:        x1   beta   0.8   0.7841 0.0462 -0.0159       -1.982   -0.3436
    #> 2:        x2   beta  -0.6  -0.5836 0.0413  0.0164       -2.726    0.3963
    #> 3:     ASC_2    asc   0.5   0.4730 0.0435 -0.0270       -5.409   -0.6222
    #> 4:     ASC_3    asc  -0.5  -0.5379 0.1082 -0.0379        7.571   -0.3498
    #> 5:  Sigma_11  sigma   1.0   1.0000 0.0000  0.0000        0.000        NA
    #> 6:  Sigma_21  sigma   0.5   0.4109 0.1341 -0.0891      -17.817   -0.6642
    #> 7:  Sigma_22  sigma   1.5   1.4108 0.2770 -0.0892       -5.949   -0.3221
    #>    lower_ci upper_ci covers
    #>       <num>    <num> <lgcl>
    #> 1:   0.6937   0.8746   TRUE
    #> 2:  -0.6645  -0.5028   TRUE
    #> 3:   0.3878   0.5581   TRUE
    #> 4:  -0.7499  -0.3258   TRUE
    #> 5:   1.0000   1.0000   TRUE
    #> 6:   0.1480   0.6738   TRUE
    #> 7:   0.8678   1.9537   TRUE

In this simulated example the posterior summaries line up with the parameters that generated the data. As with the frequentist recovery vignettes, this is an illustration rather than a proof: what matters in repeated use is posterior coverage, mixing, and sensitivity to prior and normalization choices.

## A quick MCMC diagnostic

Posterior draws are stored on the fit, so the usual diagnostics are a line away. A trace plot of the price coefficient should look like a stationary “fuzzy caterpillar”:

    beta_draws <- fit$draws$beta
    plot(beta_draws[, "x2"], type = "l", col = "steelblue",
         xlab = "iteration", ylab = expression(beta[x2]),
         main = "Posterior trace: price coefficient")
    abline(h = sim$true_params$beta[2], col = "red", lwd = 2)

The red line marks the true value; the chain should hover around it. The same diagnostics used for the hierarchical models apply to any draws matrix: `rhat()` computes the split R-hat (values near 1 are consistent with convergence) and `ess()` the rank-normalized bulk and tail effective sample sizes — the number of effectively independent draws behind each posterior summary:

    rhat(fit$draws$beta)
    #>     x1     x2  ASC_2  ASC_3
    #> 1.0043 1.0198 0.9994 1.0179
    ess(fit$draws$beta)
    #>         bulk   tail
    #> x1    152.22 516.25
    #> x2    124.65 447.31
    #> ASC_2  88.24 204.13
    #> ASC_3  55.35  67.19
    mcse(fit$draws$beta)
    #>       x1       x2    ASC_2    ASC_3
    #> 0.003741 0.003696 0.004627 0.014545

On a single chain these are necessary rather than sufficient checks; the stronger practical check in the current API is to repeat the complete fit with several independent seeds, then compare posterior summaries and pass equal-length draw matrices as a list to the same diagnostic functions, for example `rhat(list(fit$draws$beta, fit_seed_2$draws$beta), rank = TRUE)`. Because the sampler’s initial state is internal rather than user-controlled, this is not equivalent to deliberately overdispersed starting values. Longer runs, posterior stability across seeds, and prior sensitivity should therefore be reported together.

## What this fit does not provide

`choicer_mnp` is a posterior-draws object, not a `choicer_fit`. It provides coefficient and covariance summaries, recovery tools, and access to the raw and identified draws. It does not currently implement `predict()`, elasticities, diversion ratios, logsum welfare, `logLik()`, AIC, or BIC. General MNP counterfactual prediction requires multivariate-normal probability evaluation and an explicit rule for transporting the covariance structure to a new choice set; choicer does not silently impose either.

For a reportable MNP application, document the base alternative and per-draw scale normalization, the prior on the raw covariance and its implications after normalization, the alternative-specific variation that identifies covariance, trace/ESS/MCSE evidence, stability across independent seeds, and sensitivity to the prior and base alternative. Covariance conclusions are often the first to move when the identifying variation is weak.

The full derivations — the data-augmented Gibbs sampler, the McCulloch-Rossi normalization, and the identification discussion — are in the math companion. For a faster, frequentist alternative with the same post-estimation toolkit, see the multinomial logit and mixed logit vignettes. For the hierarchical extension — respondent-level random tastes and partially-pooled alternative effects, in both logit and probit flavors — see the hierarchical Bayes vignette.

## References

Albert, J. H. and Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association, 88(422), 669-679.

Keane, M. P. (1992). A note on identification in the multinomial probit model. Journal of Business & Economic Statistics, 10(2), 193-200.

McCulloch, R. and Rossi, P. E. (1994). An exact likelihood analysis of the multinomial probit model. Journal of Econometrics, 64(1-2), 207-240.

Rossi, P. E., Allenby, G. M. and McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley.

Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press, Chapter 5.