Model Fitting

Rprobit implements a number of different estimation procedures based on the minimization of a number of different criterion functions:

  1. The full probit negative log likelihood
  2. Different versions of composite marginal likelihood (CML) based on choice probabilities for pairs of choices.

Here the likelihood for the multinomial probit (MNP) model is obtained in the following way: It is assumed that the observations for different individuals are independent. Thus the log likelihood is the sum of the log likelihood contributions of the individual deciders. For each decider the choice taken in the tt-th decision is given by the random utility formulation as

Untj=Vntj+entj,n=1,...,N,t=1,...,T,j=1,...,J,Vntj=Xntj′β+Zntj′γ,γ∼N(0,Ω),ent=(ent1⋮entJ)∼N(0,Σ)\begin{array}{lll} U_{ntj} & = & V_{ntj} + e_{ntj}, n = 1,...,N, t=1,...,T, j=1,...,J, \\ V_{ntj} & = & X_{ntj}' \beta + Z_{ntj}'\gamma, \gamma \sim {\mathcal N}(0, \Omega), \\ e_{nt} & =& \begin{pmatrix} e_{nt1} \\ \vdots \\ e_{ntJ} \end{pmatrix} \sim {\mathcal N}(0,\Sigma) \end{array}

Details on the formulation of the model can be found in the vignette Model Definition.

Given this model the probability that decider nn in her tt-th choice occasion chooses alternative jj equals the probability

P(ynt=j)=P(Untk≤Untj∀k)=P(Untk−Untj≤0∀k)=P(entk−entj≤Vntj−Vntk∀k).{\mathbb P}(y_{nt} = j) = {\mathbb P}(U_{ntk} \le U_{ntj} \forall k) = {\mathbb P}( U_{ntk}-U_{ntj} \le 0 \forall k) = {\mathbb P}( e_{ntk}-e_{ntj} \le V_{ntj}-V_{ntk} \forall k). Since entke_{ntk} are multivariate normally distributed, the choice probability is obtained by evaluating the Gaussian CDF.

In a panel situation, where more than one choice occasion per decider is recorded, the different choices may be conditionally (on the regressors) independent. If Ω>0\Omega>0 such that there are random coefficients, the choices are linked. Consequently the choice occasions need to be modeled jointly. The probability of the joint vector of TT choices then obtains as

P(ynt=j,t=1,...,T)=P(entk−entj≤Vntj−Vntk,k=1,...,J,t=1,...,T,k≠j).{\mathbb P}(y_{nt} = j, t=1,...,T) = {\mathbb P}( e_{ntk}-e_{ntj} \le V_{ntj}-V_{ntk} , k=1,...,J , t=1,...,T, k\ne j).

For TT choice occasions with JJ choices each then a T(J−1)T(J-1) dimensional Gaussian CDF needs to be evaluated. This evaluation is numerically costly.

Therefore analytic approximations of the Gaussian CDF, which is used extensively in all criterion functions, as well as the corresponding gradients, are used. The package implements the Solow-Joe (SJ, Solow (1990)) approach. The approximations reduce the numerical load as compared to the usual evaluation of the Gaussian CDF. For results on the properties of the estimates using different approximations see (Batram and Bauer 2019).

The full probit likelihood evaluates the choice probabilities for all choice occasions of a decider jointly. If no mixed coefficients are present then the joint probability is the product of the marginal probabilities and the CML approaches and the full likelihood are identical (for appropriate weighting in the case of unbalanced panels).

Randomly mixed coefficients, however, introduce correlation between the choices such that the evaluation of the joint probability requires an evaluation of a potentially high dimensional normally distributed random variable. Numerically this is costly.

The CML alleviates these problems by replacing the full likelihood with easier to calculate criterion functions, see (Varin and Vidoni 2005, 2006). The most often used form is the so called pairwise likelihood which is based on representing pairs of choices in the criterion function. When picking one pair of observations from each decider leads to a criterion function with the usual properties for the reduced data set. As this is so for every choice of the pair, the same holds for a sum of such pairwise likelihoods.

Types of CML

Using pairwise CMLs leads to a criterion function that trades off numerical load with statistical efficiency. The pairwise probabilities are faster to calculate, but there is no guarantee that the CML based approach achieves the Cramer Rao lower bound. In the literature it is shown using simulations that in some cases the price in terms of asymptotic variance is small while the computational gain is large (cf. Bhat and Sidharthan (2011)).

Depending on which pairs are chosen, the tradeoff between numerical speed and variance is different. The standard choices are:

  1. cml_type = "FP": full pairwise, includes all pairs of choices.
  2. cml_type = "AC": adjacent pairwise chained. This uses all adjacent pairs plus the pair of the first and the last choice. This structure is aimed at cyclical situations such as surveys over all weekdays, such that Sunday and Monday are adjacent days.
  3. cml_type = "AP": adjacent pairwise.

Optimization Parameters

The optimization uses the R function nlm for the numerical search. This function contains many input arguments that can be specified also for the fit function of Rprobit. Amongst the most prominent of those is the argument print.level which controls the output presented in the screen. The most often used arguments are print.level = 0 (no output) and print.level = 2 (output per iteration).

Initial Estimates

The numerical search needs to be started at an initial estimate. A numeric parameter vector can be supplied through theta_init, which is useful for simulated data sets where the true underlying parameter vector is known. If theta_init = NULL, Rprobit applies its initialization procedure.

In Rprobit a strategy is used, that is described in (Oelschläger 2026): The fast initialization method uses the structure of the problem to obtain a rough guess, which can be shown to be consistent for the true parameter vector under appropriate assumptions. Simulations in (Oelschläger 2026) document that the approach leads to numerical gains compared to random initialization.

Example

Conversion of the Train data set (Croissant 2020) into a Rprobit object is described in the vignette Data and Simulation. The model can be set up as described in detail in Modeling Setup.

data("Train", package = "mlogit")

train_model <- setup(
  formula = choice ~ price + time + change | 0 | comfort,
  df = Train,
  re = c("time"),
  column_decider_id = "id",
  column_occasion_id = "choiceid"
)
train_model$mod$H_b <- train_model$mod$H_b[, -1]
train_model$mod$f_b[1] <- -0.001
train_model$mod$H_L <- matrix(c(0, 0, 1), ncol = 1)
train_model$mod$f_L <- c(0, 0, 0)

print(train_model)
#> <Rprobit object>
#> 
#> Deciders: 235
#> Choice occasions: from 5 to 19
#> Formula: choice ~ price + time + change | 0 | comfort
#> Random effects: time
#> 
#> Utility equations:
#> 
#>   U_A = b[1] * price_A + b[2] * change_A + b[3] * comfort_A + b[5] * time_A + eps[1]
#>   U_B = b[1] * price_B + b[2] * change_B + b[4] * comfort_B + b[5] * time_B + eps[2]

Having set up the model we may estimate it using the method fit for the R6 object:

start_time <- Sys.time()
train_model$fit(cml_type = "probit")
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.
end_time <- Sys.time()
elapsed_time <- end_time - start_time
elapsed_time
#> Time difference of 40.47827 secs

The initialization procedure finds initial values that are already close to the optimum:

rbind(
  "init" = train_model$theta_init,
  "fit" = train_model$theta_fit
) |> round(2)
#>       [,1]  [,2]  [,3]  [,4] [,5] [,6]
#> init -0.23 -0.68 -0.70 -0.02 0.09 1.39
#> fit  -0.23 -0.66 -0.68 -0.02 0.03 1.02

The properties of the estimated model are stored in the object and can be printed:

summary(train_model)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ price + time + change | 0 | comfort 
#> 
#> Data:
#>   Deciders: 235 
#>   Choice occasions: 2929 
#>   Choice occasions per decider: 5 to 19 
#>   Alternatives: 2 (A, B)
#> 
#> Response residuals:
#>                A         B
#> Min    -0.984846 -0.976685
#> 1Q     -0.398610 -0.407042
#> Median  0.008985 -0.008985
#> 3Q      0.407042  0.398610
#> Max     0.976685  0.984846
#> 
#> Coefficients:
#>                    Estimate Std. Error z value Pr(>|z|)    
#> price            -0.0010000         NA      NA       NA    
#> change           -0.2268586  0.0307086  -7.387 1.50e-13 ***
#> comfort_A        -0.6588193  0.0285801 -23.052  < 2e-16 ***
#> comfort_B        -0.6783702  0.0284627 -23.834  < 2e-16 ***
#> time             -0.0214950  0.0028834  -7.455 9.01e-14 ***
#> Omega[time,time]  0.0007551  0.0001308   5.772 7.84e-09 ***
#> Sigma[A,A]        0.0000000         NA      NA       NA    
#> Sigma[B,A]        0.0000000         NA      NA       NA    
#> Sigma[B,B]        1.0429832  0.0737552  14.141  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Random-effect covariance matrix (Omega):
#>           time
#> time 0.0007551
#> 
#> Error covariance matrix (Sigma):
#>   A     B
#> A 0 0.000
#> B 0 1.043
#> 
#> Prediction summary:
#>                                 A      B
#> Observed choice frequency  0.5032 0.4968
#> Mean predicted probability 0.5013 0.4987
#> 
#> Fit statistics:
#>   Log-likelihood: -1696 
#>   AIC:  3404 
#>   BIC:  3440
coef(train_model)
#>       price      change   comfort_A   comfort_B        time 
#> -0.00100000 -0.22685856 -0.65881933 -0.67837017 -0.02149496
confint(train_model)
#>                 2.5 %      97.5 %
#> price              NA          NA
#> change    -0.28704630 -0.16667083
#> comfort_A -0.71483534 -0.60280331
#> comfort_B -0.73415595 -0.62258440
#> time      -0.02714636 -0.01584357
vcov(train_model)
#>           price       change    comfort_A    comfort_B         time
#> price         0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> change        0 9.430177e-04 1.439521e-04 5.714975e-05 2.133815e-05
#> comfort_A     0 1.439521e-04 8.168235e-04 4.059792e-04 4.524468e-06
#> comfort_B     0 5.714975e-05 4.059792e-04 8.101225e-04 1.182383e-06
#> time          0 2.133815e-05 4.524468e-06 1.182383e-06 8.314096e-06

A different criterion function can be used for estimation:

start_time <- Sys.time()
train_2 <- fit(train_model, theta_init = NULL, print.level = 0, cml_type = "AP")
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.
end_time <- Sys.time()
elapsed_time2 <- end_time - start_time
elapsed_time2
#> Time difference of 16.60883 secs
summary(train_2)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ price + time + change | 0 | comfort 
#> 
#> Data:
#>   Deciders: 235 
#>   Choice occasions: 2929 
#>   Choice occasions per decider: 5 to 19 
#>   Alternatives: 2 (A, B)
#> 
#> Response residuals:
#>               A        B
#> Min    -0.98002 -0.97170
#> 1Q     -0.40360 -0.41210
#> Median  0.01396 -0.01396
#> 3Q      0.41210  0.40360
#> Max     0.97170  0.98002
#> 
#> Coefficients:
#>                    Estimate Std. Error z value Pr(>|z|)    
#> price            -0.0010000         NA      NA       NA    
#> change           -0.2052902  0.0474599  -4.326 1.52e-05 ***
#> comfort_A        -0.6554253  0.0543279 -12.064  < 2e-16 ***
#> comfort_B        -0.6730138  0.0561024 -11.996  < 2e-16 ***
#> time             -0.0188856  0.0024844  -7.602 2.92e-14 ***
#> Omega[time,time]  0.0010144  0.0002494   4.067 4.75e-05 ***
#> Sigma[A,A]        0.0000000         NA      NA       NA    
#> Sigma[B,A]        0.0000000         NA      NA       NA    
#> Sigma[B,B]        1.1565400  0.1414721   8.175 2.96e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Random-effect covariance matrix (Omega):
#>          time
#> time 0.001014
#> 
#> Error covariance matrix (Sigma):
#>   A     B
#> A 0 0.000
#> B 0 1.157
#> 
#> Prediction summary:
#>                                 A      B
#> Observed choice frequency  0.5032 0.4968
#> Mean predicted probability 0.5009 0.4991
#> 
#> Fit statistics:
#>   Log-likelihood: -1750 
#>   AIC:  3512 
#>   BIC:  3548

The AP criterion is much faster in this example. Because the two fits optimize different criterion functions, their criterion values are not directly comparable.

c(logLik(train_model), logLik(train_2))
#> [1] -1695.945 -1750.118

Their estimated standard errors can be compared as follows:

V_probit <- vcov(train_model)
V_ap <- vcov(train_2)

cbind(sqrt(diag(V_probit)), sqrt(diag(V_ap)))
#>                  [,1]        [,2]
#> price     0.000000000 0.000000000
#> change    0.030708593 0.047459882
#> comfort_A 0.028580124 0.054327931
#> comfort_B 0.028462651 0.056102424
#> time      0.002883417 0.002484365

Model Diagnostics

Beside the fitted model also the predictions of the choice probabilities as well as residuals are of interest. These can be obtained as follows:

res <- residuals(train_model)
apply(res, 2, mean)
#>            A            B 
#>  0.001922934 -0.001922934
probs <- predict(train_model)
plot(probs[, 1], res[, 1], xlab = "Choice probabilities", ylab = "Residuals")

The plot shows that the residuals are given as the choices minus the choice probabilities.

Additionally predicted choice probabilities can also be obtained for a new data set:

prediction_df <- Train[1:6, ]
prediction_df[, "price_A"] <- 3000
prediction_df[, "time_A"] <- 120
prediction_df[, "choice"] <- NA
prediction_data <- choice_data$new(
  df = prediction_df,
  column_decider_id = "id"
)
new_probabilities <- predict(
  train_model,
  newdata = prediction_data,
  type = "probabilities"
)
new_probabilities
#>           A         B
#> 1 0.8976261 0.1023739
#> 2 0.6594053 0.3405947
#> 3 0.5897089 0.4102911
#> 4 0.5563572 0.4436428
#> 5 0.5563572 0.4436428
#> 6 0.3578979 0.6421021

References

Batram, M., and D. Bauer. 2019. “On Consistency of the MACML Approach to Discrete Choice Modelling.” Journal of Choice Modelling 30: 1–16. https://doi.org/10.1016/j.jocm.2018.10.001.
Bhat, C., and R. Sidharthan. 2011. “A Simulation Evaluation of the Maximum Approximate Composite Marginal Likelihood (MACML) Estimator for Mixed Multinomial Probit Models.” Transportation Research Part B: Methodological 45 (7): 940–53.
Croissant, Y. 2020. “Estimation of Random Utility Models in R: The mlogit Package.” Journal of Statistical Software 95 (11). https://doi.org/10.18637/jss.v095.i11.
Oelschläger, L. 2026. “Overcoming Challenges in Modeling Choice Behavior Heterogeneity.” PhD thesis, Bielefeld University. https://pub.uni-bielefeld.de/record/3014719.
Solow, A. 1990. “A Method for Approximating Multivariate Normal Orthant Probabilities.” Journal of Statistical Computation and Simulation 37 (3-4): 225–29.
Varin, C., and P. Vidoni. 2005. “A Note on Composite Likelihood Inference and Model Selection.” Biometrika 92 (3): 519–28.
Varin, C., and P. Vidoni. 2006. “Pairwise Likelihood Inference for Ordinal Categorical Time Series.” Computational Statistics & Data Analysis 51 (4): 2365–73.