Model definition

The probit model

The probit model¹ is a regression-type model where the dependent variable only takes a finite number of values and the error term is normally distributed (Agresti 2015). Its purpose is to estimate the probability that the dependent variable takes a certain, discrete value. The most common application are discrete choice scenarios. The dependent variable here is one of finitely many and mutually exclusive alternatives, and explanatory variables typically are characteristics of the deciders or the alternatives.

To be concrete, assume that we possess data of N decision makers which choose between J ≥ 2 alternatives² at each of T choice occasions³. Specific to each decision maker, alternative and choice occasion, we furthermore observe P choice attributes that we use to explain the choices. The continuous choice attributes cannot be linked directly to the discrete choices but must take a detour over a latent variable. In the discrete choice setting, this variable can be interpreted as the decider’s utility of a certain alternative. Decider n’s utility U_ntj for alternative j at choice occasion t is modeled as

for n = 1, …, N, t = 1, …, T and j = 1, …, J, where

• X_ntj is a (column) vector of P characteristics of j as faced by n at t,

• β ∈ ℝ^P is a vector of coefficients,

• and (ϵ_nt:) = (ϵ_nt1, …, ϵ_ntJ)′ ∼ MVN_J(0, Σ) is the model’s error term vector for n at t, which in the probit model is assumed to be multivariate normally distributed with zero mean and covariance matrix Σ.

Now let y_nt = j denote the event that decision maker n chooses alternative j at choice occasion t. Assuming utility maximizing behavior of the decision makers⁴, the decisions are linked to the utilities via

In the ordered probit case, the concept of decider’s having separate utilities for each alternative is no longer natural (Train 2009). Instead, we model only a single utility value per decider n and choice occasion t, which we interpret as the “level of association” that n has with the choice question. The utility value falls into discrete categories, which in turn are linked to the ordered alternatives j = 1, …, J. Formally, with end points γ₀ = −∞ and γ_J = +∞, and thresholds (γ_j)_{j = 1, …, J − 1}. To ensure monotonicity of the thresholds, we rather estimate logarithmic threshold increments d_j with γ_j = ∑_{i = 1, …, j}exp d_i, j = 1, …, J − 1.

Choice behavior heterogeneity

Note that the coefficient vector β is constant across decision makers. This assumption is too restrictive for many applications.⁵ Heterogeneity in choice behavior can be modeled by imposing a distribution on β such that each decider can have their own preferences.

Formally, we define β = (α, β_n), where α are P_f coefficients that are constant across deciders and β_n are P_r decider-specific coefficients. Consequently, P = P_f + P_r. Now if P_r > 0, β_n is distributed according to some P_r-variate distribution, the so-called mixing distribution.

Choosing an appropriate mixing distribution is a notoriously difficult task of the model specification. In many applications, different types of standard parametric distributions (including the normal, log-normal, uniform and tent distribution) are tried in conjunction with a likelihood value-based model selection, cf., Train (2009), Chapter 6. Instead, {RprobitB} implements the approach of (Oelschläger and Bauer 2020) to approximate any underlying mixing distribution by a mixture of (multivariate) Gaussian densities. More precisely, the underlying mixing distribution g_Pr for the random coefficients (β_n)_n is approximated by a mixture of P_r-variate normal densities ϕ_Pr with mean vectors b = (b_c)_c and covariance matrices Ω = (Ω_c)_c using C components, i.e.

Here, (s_c)_c are weights satisfying 0 < s_c ≤ 1 for c = 1, …, C and ∑_cs_c = 1. One interpretation of the latent class model is obtained by introducing variables z = (z_n)_n, allocating each decision maker n to class c with probability s_c, i.e.

We call the resulting model the latent class mixed multinomial probit model. Note that the model collapses to the (normally) mixed multinomial probit model if P_r > 0 and C = 1, to the multinomial probit model if P_r = 0 and to the binary probit model if additionally J = 2.

Model normalization

As is well known, any utility model needs to be normalized with respect to level and scale in order to be identified (Train 2009). Therefore, we consider the transformed model

n = 1, …, N, t = 1, …, T and j = 1, …, J − 1, where (choosing J as the reference alternative) Ũ_ntj = U_ntj − U_ntJ, X̃_ntj = X_ntj − X_ntJ, and ϵ̃_ntj = ϵ_ntj − ϵ_ntJ, where (ϵ̃_nt:) = (ϵ̃_nt1, ..., ϵ̃_{nt(J − 1)})′ ∼ MVN_J − 1(0, Σ̃) and Σ̃ denotes a covariance matrix with the top-left element restricted to one.⁶

Parameter labels

In {RprobitB}, the probit model parameters are saved as an RprobitB_parameter object. Their labels are consistent with their definition in this vignette. For example:

RprobitB:::RprobitB_parameter(
  P_f = 1,
  P_r = 2,
  J = 3,
  N = 10,
  C = 2, # the number of latent classes
  alpha = c(1), # the fixed coefficient vector of length 'P_f'
  s = c(0.6, 0.4), # the vector of class weights of length 'C'
  b = matrix(c(-1, 1, 1, 2), nrow = 2, ncol = 2),
  # the matrix of class means as columns of dimension 'P_r' x 'C'
  Omega = matrix(c(diag(2), 0.1 * diag(2)), nrow = 4, ncol = 2),
  # the matrix of class covariance matrices as columns of dimension 'P_r^2' x 'C'
  Sigma = diag(2), # the differenced error term covariance matrix of dimension '(J-1)' x '(J-1)'
  # the undifferenced error term covariance matrix is labeled 'Sigma_full'
  z = rep(1:2, 5) # the vector of the allocation variables of length 'N'
)
#> alpha : 1
#> 
#> C : 2
#> 
#> s : double vector of length 2 
#> 0.6 0.4
#> 
#> b : 2 x 2 matrix of doubles 
#>      [,1] [,2]
#> [1,]   -1    1
#> [2,]    1    2
#> 
#> 
#> Omega : 4 x 2 matrix of doubles 
#>      [,1] [,2]
#> [1,]    1  0.1
#> [2,]    0    0
#> [3,]    0    0
#> [4,]    1  0.1
#> 
#> 
#> Sigma : 2 x 2 matrix of doubles 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> 
#> Sigma_full : 3 x 3 matrix of doubles 
#>      [,1] [,2] [,3]
#> [1,]    2    1    1
#> [2,]    1    2    1
#> [3,]    1    1    1
#> 
#> 
#> beta : 2 x 10 matrix of doubles 
#>      [,1] [,2]  [,3] ... [,10]
#> [1,] 1.14 0.82 -1.14 ...  0.34
#> [2,] 1.11 2.03  0.37 ...  2.21
#> 
#> 
#> z : double vector of length 10 
#> 1 2 1 ... 2
#> 
#> d : NA

Mind that the matrix Sigma_full is not unique and can be any matrix that results into Sigma after the differencing, see the non-exported function RprobitB:::undiff_Sigma().