Model Definition

The Multinomial Probit Model

The multinomial probit model¹ provides a mathematical representation of process of choosing one out of a finite number of mutually exclusive alternatives depending on a number of regressor variables characterizing the alternatives and the decider. It possesses the interpretation of a random utility maximization (RUM) model, whereby a decider assigns a utility value to each alternative and picks the one with highest utility. As the modeler does not have full information, the model represents the utility using a systematic part, depending on regressor variables, plus an additive error term, which for probit models is normally distributed (Agresti 2015).

The model output are conditional choice probabilities depending on the regressor values.² The regressor variables potentially consist of characteristics describing the alternatives and the deciders. To be concrete, assume that we possess data of $N$ decision makers which choose between $J \geq 2$ alternatives at each of $T$ choice occasions³. One example can be the choice of one out of a number of cars, both conventional and electric, for the next purchase. Specific to each decision maker, alternative and choice occasion, we observe $P$ choice attributes (like the purchase price, for example) that we use to explain the choices. Decider $n$’s utility $U_{ntj}$ for alternative $j$ at choice occasion $t$ is modeled as

$$U_{ntj} = X_{ntj}'\beta + S_{nt}' \gamma_j + Z_{ntj}'\delta_j + \alpha_j + \epsilon_{ntj}$$

for $n=1,\dots,N$, $t=1,\dots,T$ and $j=1,\dots,J$, where

• $X_{ntj}$ ($Z_{ntj}$) is a (column) vector of $P$ ($P_Z$) characteristics of alternative $j$ as faced by $n$ at choice occasion $t$,

• $S_{nt}$ is a (column) vector of $P_S$ characteristics of the decider $n$ at occasion $t$,

• $\beta \in {\mathbb R}^{P}, \gamma_j \in {\mathbb R}^{P_S}, \delta_j \in {\mathbb R}^{P_Z}$ is a vector of coefficients,

• $\alpha_j \in {\mathbb R}$ are called alternative specific constants (ASCs),

• and $(\epsilon_{nt:}) = (\epsilon_{nt1},\dots,\epsilon_{ntJ})' \sim \text{MVN}_{J} (0,\Sigma)$ 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 $\Sigma$. Details on how such data sets are represented in Rprobit are given in the vignette Data and Simulation.

The three types of regressor variables are called

• type I variables ($X_{ntj}$) that vary over alternatives but influence the utility with a common coefficient (that is, not depending on the alternative),

• type II variables ($S_{nt}$) that do not vary over alternatives and influence the utility with coefficients differing over alternatives,

• type III variables ($Z_{ntj}$) that vary over alternatives but influence the utility with coefficients differing over alternatives.

Now let $y_{nt}=j$ denote the event that decision maker $n$ chooses alternative $j$ at her $t$-th choice occasion. Assuming utility maximizing behavior of the decision makers, the decisions are linked to the utilities via⁴

$$y_{nt} = {\arg \max}_{j = 1,\dots,J} U_{ntj}.$$

Model Formula

In Rprobit the model is represented using a formula as in the mlogit R package (Croissant 2020). The following formula describes the dependent variable choice as a function of the type I variable V1, the type II variable V2 and the type III variable V3:

formula <- choice ~ V1 | V2 | V3

Hereby the following rules apply:

• choice must be categorical variable listing the choices taken.
• Type III variables may be omitted.
• Type II variables may be omitted, if type III are also omitted. Otherwise include a 0
• More than one variable can be included in each of the three sections using addition signs such as V1 + V2 + V3.
• By default ASCs are included as type II variables. They can be excluded by adding + 0 such as V2 + 0.

Example: Choice of a New Car

Consider a household deciding between different power trains for the car to be purchased. The options are a battery powered electric vehicle (BEV) and an internal combustion engine (ICE).

The alternative specific variables may include:

• price price
• driving range range
• horse power equivalent PS
• expected maintenance price per year Main

All of these depend on the alternative at the time of buying. The driving range coefficient for the BEV will be much higher than for the ICEV, as the refilling of the ICEV is simpler.

The decider characteristics may include:

• household income HI
• household savings HS
• location (urban?) loc
• charging facility potential CF
• attitude towards environmental issues Attitude

The corresponding formula might be assuming that Attitude is not measured:

formula <- choice ~  PS + Main | HI + HS + loc + CF + 0 | range

The tradeoff between the purchase price and the driving range of a car may depend on the commuting habits of the deciders that are not represented well in the data.

Identifiability

When maximizing utilities the level is not of importance: Adding a constant to all utilities does not change the maximizing alternative. Also the scale is not of importance: Multiplying all utilities with a positive constant leaves the maximizing alternative indifferent. This implies identifiability issues. (Train 2009)

For type I variables the difference between two utilities $U_{ntj}-U_{ntk}$ depends on the difference $X_{ntj}-X_{ntk}$ between the regressor variables. Therefore variables that do not vary over alternatives (resulting in $X_{ntj}-X_{ntk}=0$) cannot be used as type I variables. Additionally, $\beta$ needs to be identified from $X_{ntj}'\beta$. Thus the type I variables must not be multicollinear.

For each type II variable the model includes $J$ coefficients $\gamma_j, j=1,...,J$. Subtracting $S_{nt}'\gamma_1$ from all utilities eliminates the coefficients corresponding to the first alternative from the model. The remaining $(J-1) \cdot P_S$ coefficients for the type II variables are identifiable, if there are no multicollinearities between the type II variables. Here the choice of the first variable is arbitrary and any base alternative (=reference) can be chosen.

For type III variables we obtain $Z_{ntj}'\delta_j - Z_{ntk}'\delta_k$ in the utility differences. It follows that here all $J \cdot P_Z$ coefficients are identified, if $Z_{ntj}$ vary over alternatives and are not multicollinear.

The same issue arises for the error terms: Subtracting the error $\epsilon_{nt1}$ from all utilities does not change the maximizing alternative. Thus we can without restriction of generality assume that $\epsilon_{nt1} = 0$ (or pick any alternative as a reference alternative in this respect).

Finally for the identification of the scale one may either restrict one coefficient to unity (or minus unity). Or one might restrict one variance to a particular value (typically 1 or 0.5). For details see Train (2009) Section 5.2. Numerically fixing one coefficient to $\pm 1$ is more stable than fixing a variance.

Model Unification

While the specification of type I, type II and type III variables is convenient for model specification and data storage, it is inconvenient in the estimation phase. Thus Rprobit uses a unified framework, recoding type II and type III variables into a set of type I variables. This is exemplified below. The choice data object (described in detail in the vignette Data and Simulation) provides the method build_design to perform the recoding according to a model formulation.

Choice Behavior Heterogeneity

Note that all coefficient vectors $\beta, \gamma_j, \delta_j$ above are constant across decision makers. This assumption is too restrictive for many applications.⁵ Unobserved heterogeneity in choice behavior can be modeled by imposing a distribution on any element of $\beta, \gamma_j, \delta_j$ such that each decider can have their own preferences.

Such models use a random coefficient setup wherein the random coefficients are assumed to be drawn independent of all regressor variables and the errors. Mathematically this implies a mixing of the implied choice probabilities (the expectation over the random terms mathematically leads to an integration). Rprobit currently only implements normally distributed random parameters, that can be correlated across coefficients.

More precisely, in the unified model representation using only type I variables some of the components of $\beta$ can be mixed using a normal distribution for mixing. We assume that the mixed coefficients appear in the trailing part of $\beta$, that is $\beta = [\beta_f',\beta_r'], \beta_r \in {\mathbb R}^R, \beta_r \sim {\mathcal N}(b_r, \Omega), \beta_f = b_f \in {\mathbb R}^{K-R}$.

It is assumed that all coefficients corresponding to one original type II or type III variables are either mixed or not mixed. Thus the mixed coefficients can be specified by providing the names of their corresponding variables. For example in the model formula <- choice ~ V1 | V2 | V3 the coefficients to variable V1 and V3 are defined as mixing via

re <- c("V1", "V3")

Model Parameterization

According to the unification described above we achieve model specification using only type I variables as

$$U_{ntj} = X_{ntj}'\beta + \epsilon_{ntj}, j = 1,...,J, X_{ntj} , (\epsilon_{nt:}) \in {\mathbb R}^J, \beta \in {\mathbb R}^K.$$

Some of the components of $\beta$ can be mixed using a normal distribution for mixing. We assume that the mixed coefficients appear in the trailing part of $\beta$, that is $\beta = [\beta_f',\beta_r'], \beta_r \in {\mathbb R}^R, \beta_r \sim {\mathcal N}(b_r, \Omega), \beta_f = b_f \in {\mathbb R}^{K-R}$.

Therefore the parameterization of the model requires parameters for the vector $b = (b_f',b_r')' \in {\mathbb R}^K$ and the two variance matrices $\Omega \in {\mathbb R}^{R \times R}$ and $\Sigma \in {\mathbb R}^{J \times J}$. The parameterization has to conform to the restriction that variance matrices are positive semidefinite. Additionally the parameterization must ascertain identifiability of the parameters.

The parameter vector $b \in {\mathbb R}^{K}$ is parameterized using a linear affine mapping, to allow for the imposition of some restrictions:

$$b = H_b \theta_b + f_b.$$

This allows to incorporate linear equality constraints on the parameters. For instance, restricting $b_1 = 1, b_2 = 0, b_3 = - b_4$ with $b_5$ free is achieved using

$$H_b = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ -1 & 0 \\ 0 & 1 \end{pmatrix}, f_b = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}.$$

A common device to ensure positive semidefiniteness of matrices is to use Cholesky factors. For each positive semidefinite matrix $\Sigma \in {\mathbb R}^{J \times J}$ there exists a lower triangular matrix $L \in {\mathbb R}^{J \times J}$ such that $L L' = \Sigma$.

We chose to parameterize the vectorization of the lower half of the Cholesky matrices using a linear affine mapping:

$$\operatorname{vech}(L) = H_L \theta_L + f_L$$

where $\operatorname{vech}$ denotes columnwise vectorization of the lower triangular matrix (including the diagonal). Thus for $J=3$, for example, the lower triangular Cholesky factor includes 6 entries:

$$\begin{align*} L = \begin{pmatrix} \theta_1 & 0 & 0 \\ \theta_2 & \theta_4 & 0 \\ \theta_3 & \theta_5 & \theta_6 \end{pmatrix}. \end{align*}$$

Here $f_L$ can be used to fix entries for identifiability reasons. For example, for $J=3$ parameterizing a diagonal matrix where the (1,1) entry is fixed to equal $\Sigma_{1,1}=1$ can be done using

$$\begin{align*} H_L &= \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{pmatrix}, f_L = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}. \end{align*}$$

Thus there are two parameters $\theta_1$ and $\theta_2$. We do not impose non-negativity of these parameters, which leads to a non-identifiability as without this restriction on the diagonal entries the Cholesky factors are not unique. The tradeoff made here is the inclusion of multiple well separated global optima in favor of being able to leave out parameter non-negativity restrictions.

Model Object

Rprobit collects all information related to the model in a mod object. This object is implemented as a R6 class. See here for the perils of R6 classes.

The following example parameterizes a model for mode choice depending on the price of the alternative modes (used as a type I variable) and time (used as a type III variable) featuring $J = 3$ alternatives. The ASC are mixed using a normal distribution:

formula <- choice ~ price | 1 | time
re <- "ASC"
m_basic <- mod$new(formula, re, alt = 3)

Printing the object provides an overview of the specification:

print(m_basic)
#> <mod object>
#> 
#> Formula: choice ~ price | 1 | time
#> Random effects: ASC
#> Alternative names: A, B, C
#> Coefficients: 4 fixed, 2 random
#> Free parameters total: 10
#> 
#> Utility equations:
#> 
#>   U_A = b[1] * price_A + b[2] * time_A + eps[1]
#>   U_B = b[1] * price_B + b[3] * time_B + b[5] + eps[2]
#>   U_C = b[1] * price_C + b[4] * time_C + b[6] + eps[3]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b
#>   length(theta_b): 6
#>   H_b (6 x 6) =
#>             theta_1 theta_2 theta_3 theta_4 theta_5 theta_6
#> b[1]:price        1       0       0       0       0       0
#> b[2]:time_A       0       1       0       0       0       0
#> b[3]:time_B       0       0       1       0       0       0
#> b[4]:time_C       0       0       0       1       0       0
#> b[5]:ASC_B        0       0       0       0       1       0
#> b[6]:ASC_C        0       0       0       0       0       1
#> 
#> 
#>   O = vech(chol(Omega)) = H_O theta_O
#>   length(theta_O): 2
#>   H_O (3 x 2) =
#>               theta_7 theta_8
#> ASC_B_:_ASC_B       1       0
#> ASC_C_:_ASC_B       0       0
#> ASC_C_:_ASC_C       0       1
#> 
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L + f_L
#>   length(theta_L): 2
#>   H_L (6 x 2) =
#>       theta_9 theta_10
#> A_:_A       0        0
#> B_:_A       0        0
#> C_:_A       0        0
#> B_:_B       0        0
#> C_:_B       1        0
#> C_:_C       0        1
#>   f_L (6 x 1) =
#>        
#> A_:_A 0
#> B_:_A 0
#> C_:_A 0
#> B_:_B 1
#> C_:_B 0
#> C_:_C 0

We can see that six coefficient parameters (length_theta_b = 6) are free, the two ASCs are random leading to a $2 \times 2$ covariance matrix $\Omega \in {\mathbb R}^{2 \times 2}$ that is modeled as a diagonal matrix using two parameters (length_theta_O = 2). Finally the variance $\Sigma$ uses two parameters (length_theta_L = 2):

$$\begin{align*} H_L &= \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}, f_L = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} \Rightarrow L = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \theta_1 & \theta_2 \end{pmatrix} \Rightarrow \Sigma = LL' = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \theta_1 & \theta_1^2 + \theta_2^2 \end{pmatrix}. \end{align*}$$

There are more options to initialize a mod object:

alt_names <- c("train", "bus", "car")

m <- mod$new(
  formula = formula, re = re, alt = 3,
  alt_names = alt_names,            # alternative labels
  ref_alt = "train",                # reference alternative (for type-2 effects)
  corr_re = FALSE,                  # structure for Omega
  error_struc = "diag"              # structure for Sigma
)

print(m)
#> <mod object>
#> 
#> Formula: choice ~ price | 1 | time
#> Random effects: ASC
#> Alternative names: train, bus, car
#> Coefficients: 4 fixed, 2 random
#> Free parameters total: 10
#> 
#> Utility equations:
#> 
#>   U_train = b[1] * price_train + b[2] * time_train + eps[1]
#>   U_bus = b[1] * price_bus + b[3] * time_bus + b[5] + eps[2]
#>   U_car = b[1] * price_car + b[4] * time_car + b[6] + eps[3]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b
#>   length(theta_b): 6
#>   H_b (6 x 6) =
#>                 theta_1 theta_2 theta_3 theta_4 theta_5 theta_6
#> b[1]:price            1       0       0       0       0       0
#> b[2]:time_train       0       1       0       0       0       0
#> b[3]:time_bus         0       0       1       0       0       0
#> b[4]:time_car         0       0       0       1       0       0
#> b[5]:ASC_bus          0       0       0       0       1       0
#> b[6]:ASC_car          0       0       0       0       0       1
#> 
#> 
#>   O = vech(chol(Omega)) = H_O theta_O
#>   length(theta_O): 2
#>   H_O (3 x 2) =
#>                   theta_7 theta_8
#> ASC_bus_:_ASC_bus       1       0
#> ASC_car_:_ASC_bus       0       0
#> ASC_car_:_ASC_car       0       1
#> 
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L + f_L
#>   length(theta_L): 2
#>   H_L (6 x 2) =
#>               theta_9 theta_10
#> train_:_train       0        0
#> bus_:_train         0        0
#> car_:_train         0        0
#> bus_:_bus           1        0
#> car_:_bus           0        0
#> car_:_car           0        1
#>   f_L (6 x 1) =
#>                
#> train_:_train 1
#> bus_:_train   0
#> car_:_train   0
#> bus_:_bus     0
#> car_:_bus     0
#> car_:_car     0

This adds labels to the alternatives and changes the error covariance matrix to a diagonal structure. For the recoding ref_alt defines the reference alternative, for which no ASC or type II regressors are generated to avoid non-identifiability. The correlation structure of $\Omega$ is set to FALSE leading to a diagonal matrix. The error variance matrix is chosen diagonal. The available error structures are:

• "full": the first Cholesky column is fixed to zero, L[2,2] is fixed to one, and the remaining entries are free
• "diag": the first Cholesky diagonal entry is fixed to one and the remaining diagonal entries are free
• "ident": $\Sigma$ is proportional to the identity matrix, with one shared free Cholesky diagonal parameter
• "fixed": $\Sigma$ is fixed to the identity matrix.

Working with the Model Objects

Having defined objects, one wants to manipulate the models. To do this, we must know the fields of the model object. These can be grouped into fields that are read-only and fields that can be adapted.

Read-only fields are:

• alt number of alternatives
• b_count vector of two integers, number of fixed and random coefficients
• var_names list containing the names of fixed and random coefficients
• length_theta number of free model parameters, together with length_theta_b, length_theta_O, and length_theta_L splitting it into the parts corresponding to $\beta, \Omega, \Sigma$.
• mod_list list containing the model representation used by the numerical routines.

These cannot be changed directly but are inferred from the adjustable fields:

• formula: model formula
• re: vector of variables with random coefficients
• alt_names: vector of strings labeling the alternatives
• ref_alt string identifying reference alternative
• H_b, f_b, H_O, f_O, H_L, f_L: matrices implementing the affine linear restrictions for the various matrices.
• suppress_warnings: Boolean.

The last field suppresses warnings issued by the object. Changing model-defining fields such as formula, re, alt_names, or ref_alt can automatically adapt dependent restriction matrices and vectors. By default, the object warns about these adaptations.

To demonstrate this behavior, consider adding a random coefficient for price:

m$suppress_warnings <- TRUE
m$re <- c("ASC", "price")
m
#> <mod object>
#> 
#> Formula: choice ~ price | 1 | time
#> Random effects: ASC, price
#> Alternative names: train, bus, car
#> Coefficients: 3 fixed, 3 random
#> Free parameters total: 11
#> 
#> Utility equations:
#> 
#>   U_train = b[1] * time_train + b[4] * price_train + eps[1]
#>   U_bus = b[2] * time_bus + b[4] * price_bus + b[5] + eps[2]
#>   U_car = b[3] * time_car + b[4] * price_car + b[6] + eps[3]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b
#>   length(theta_b): 6
#>   H_b (6 x 6) =
#>                 theta_1 theta_2 theta_3 theta_4 theta_5 theta_6
#> b[1]:time_train       1       0       0       0       0       0
#> b[2]:time_bus         0       1       0       0       0       0
#> b[3]:time_car         0       0       1       0       0       0
#> b[4]:price            0       0       0       1       0       0
#> b[5]:ASC_bus          0       0       0       0       1       0
#> b[6]:ASC_car          0       0       0       0       0       1
#> 
#> 
#>   O = vech(chol(Omega)) = H_O theta_O
#>   length(theta_O): 3
#>   H_O (6 x 3) =
#>                   theta_7 theta_8 theta_9
#> price_:_price           1       0       0
#> ASC_bus_:_price         0       0       0
#> ASC_car_:_price         0       0       0
#> ASC_bus_:_ASC_bus       0       1       0
#> ASC_car_:_ASC_bus       0       0       0
#> ASC_car_:_ASC_car       0       0       1
#> 
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L + f_L
#>   length(theta_L): 2
#>   H_L (6 x 2) =
#>               theta_10 theta_11
#> train_:_train        0        0
#> bus_:_train          0        0
#> car_:_train          0        0
#> bus_:_bus            1        0
#> car_:_bus            0        0
#> car_:_car            0        1
#>   f_L (6 x 1) =
#>                
#> train_:_train 1
#> bus_:_train   0
#> car_:_train   0
#> bus_:_bus     0
#> car_:_bus     0
#> car_:_car     0

Since suppress_warnings is set to TRUE, no warning is issued. If we set suppress_warnings to FALSE and change the model back, we get a warning:

m$suppress_warnings <- FALSE
m$re <- c("ASC")
m
#> <mod object>
#> 
#> Formula: choice ~ price | 1 | time
#> Random effects: ASC
#> Alternative names: train, bus, car
#> Coefficients: 4 fixed, 2 random
#> Free parameters total: 10
#> 
#> Utility equations:
#> 
#>   U_train = b[1] * price_train + b[2] * time_train + eps[1]
#>   U_bus = b[1] * price_bus + b[3] * time_bus + b[5] + eps[2]
#>   U_car = b[1] * price_car + b[4] * time_car + b[6] + eps[3]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b
#>   length(theta_b): 6
#>   H_b (6 x 6) =
#>                 theta_1 theta_2 theta_3 theta_4 theta_5 theta_6
#> b[1]:price            1       0       0       0       0       0
#> b[2]:time_train       0       1       0       0       0       0
#> b[3]:time_bus         0       0       1       0       0       0
#> b[4]:time_car         0       0       0       1       0       0
#> b[5]:ASC_bus          0       0       0       0       1       0
#> b[6]:ASC_car          0       0       0       0       0       1
#> 
#> 
#>   O = vech(chol(Omega)) = H_O theta_O
#>   length(theta_O): 2
#>   H_O (3 x 2) =
#>                   theta_7 theta_8
#> ASC_bus_:_ASC_bus       1       0
#> ASC_car_:_ASC_bus       0       0
#> ASC_car_:_ASC_car       0       1
#> 
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L + f_L
#>   length(theta_L): 2
#>   H_L (6 x 2) =
#>               theta_9 theta_10
#> train_:_train       0        0
#> bus_:_train         0        0
#> car_:_train         0        0
#> bus_:_bus           1        0
#> car_:_bus           0        0
#> car_:_car           0        1
#>   f_L (6 x 1) =
#>                
#> train_:_train 1
#> bus_:_train   0
#> car_:_train   0
#> bus_:_bus     0
#> car_:_bus     0
#> car_:_car     0

To impose a restriction on $\beta$, we may change the matrices H_b, f_b directly. For example, restricting $b(1)$ to equal $b[1]=1$ we may use the following commands:

m$H_b <- m$H_b[, -1]
m$f_b[1] <- 1
m
#> <mod object>
#> 
#> Formula: choice ~ price | 1 | time
#> Random effects: ASC
#> Alternative names: train, bus, car
#> Coefficients: 4 fixed, 2 random
#> Free parameters total: 9
#> 
#> Utility equations:
#> 
#>   U_train = b[1] * price_train + b[2] * time_train + eps[1]
#>   U_bus = b[1] * price_bus + b[3] * time_bus + b[5] + eps[2]
#>   U_car = b[1] * price_car + b[4] * time_car + b[6] + eps[3]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b + f_b
#>   length(theta_b): 5
#>   H_b (6 x 5) =
#>                 theta_1 theta_2 theta_3 theta_4 theta_5
#> b[1]:price            0       0       0       0       0
#> b[2]:time_train       1       0       0       0       0
#> b[3]:time_bus         0       1       0       0       0
#> b[4]:time_car         0       0       1       0       0
#> b[5]:ASC_bus          0       0       0       1       0
#> b[6]:ASC_car          0       0       0       0       1
#>   f_b (6 x 1) =
#>             
#> price      1
#> time_train 0
#> time_bus   0
#> time_car   0
#> ASC_bus    0
#> ASC_car    0
#> 
#> 
#>   O = vech(chol(Omega)) = H_O theta_O
#>   length(theta_O): 2
#>   H_O (3 x 2) =
#>                   theta_6 theta_7
#> ASC_bus_:_ASC_bus       1       0
#> ASC_car_:_ASC_bus       0       0
#> ASC_car_:_ASC_car       0       1
#> 
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L + f_L
#>   length(theta_L): 2
#>   H_L (6 x 2) =
#>               theta_8 theta_9
#> train_:_train       0       0
#> bus_:_train         0       0
#> car_:_train         0       0
#> bus_:_bus           1       0
#> car_:_bus           0       0
#> car_:_car           0       1
#>   f_L (6 x 1) =
#>                
#> train_:_train 1
#> bus_:_train   0
#> car_:_train   0
#> bus_:_bus     0
#> car_:_bus     0
#> car_:_car     0

To check that the first element is restricted to one, we can draw a random parameter vector and calculate the system corresponding to this parameter vector:

set.seed(1)
theta <- runif(m$length_theta)
(parameter <- m$build_parameter(theta))
#> $b
#>                     
#> price      1.0000000
#> time_train 0.2655087
#> time_bus   0.3721239
#> time_car   0.5728534
#> ASC_bus    0.9082078
#> ASC_car    0.2016819
#> 
#> $O
#>                            
#> ASC_bus_:_ASC_bus 0.8983897
#> ASC_car_:_ASC_bus 0.0000000
#> ASC_car_:_ASC_car 0.9446753
#> 
#> $Omega_chol
#>           ASC_bus   ASC_car
#> ASC_bus 0.8983897 0.0000000
#> ASC_car 0.0000000 0.9446753
#> 
#> $Omega
#>          ASC_bus   ASC_car
#> ASC_bus 0.807104 0.0000000
#> ASC_car 0.000000 0.8924114
#> 
#> $L
#>                        
#> train_:_train 1.0000000
#> bus_:_train   0.0000000
#> car_:_train   0.0000000
#> bus_:_bus     0.6607978
#> car_:_bus     0.0000000
#> car_:_car     0.6291140
#> 
#> $Sigma_chol
#>       train       bus      car
#> train     1 0.0000000 0.000000
#> bus       0 0.6607978 0.000000
#> car       0 0.0000000 0.629114
#> 
#> $Sigma
#>       train       bus       car
#> train     1 0.0000000 0.0000000
#> bus       0 0.4366537 0.0000000
#> car       0 0.0000000 0.3957845

Sometimes it is useful to go the other way and calculate the parameter vector corresponding to a particular system:

theta2 <- m$build_theta(
  b = parameter$b,
  Omega = parameter$Omega,
  Sigma = parameter$Sigma
)
all.equal(theta, theta2)
#> [1] TRUE

This helps for example when a given variance matrix $\Sigma$ shall be parameterized.

Restrictions on $\Omega$ and $\Sigma$ can be handled analogously. Note, that Rprobit provides much freedom with respect to choosing the parameterization of the model. The user can change the matrices and vectors in the three affine linear mappings within the dimensions implied by the model. Consequently, the user has to make sure that her parameterization is identifiable. (Train 2009) discusses in section 5.2. guidelines to achieve identification. This is supported via the method check_identifiability which performs a numeric check by calculating the column rank of the Jacobian of a normalized representation of $b$, $\Omega$, and the differenced $\Sigma$ for a randomly chosen parameter vector $\theta$.

If in the example above, for example, we change the parameterization back to estimate all entries in $\beta$ and to also include the $(2,2)$-element of $L$ as a free parameter, then the model loses its identifiability due to a lack of identifying the scale of the utility. Some problems with identifiability can be checked using the check_identifiability() function. The first model is identifiable…

m$H_b <- diag(6)
m$check_identifiability()
#> [1] TRUE

… while the second is not:

m2 <- m$clone()
m2$H_L <- diag(6)[, -c(1, 2, 3)]
m2$f_L <- matrix(0, ncol = 1, nrow = 6)
m2$check_identifiability()
#> Warning: The parameter set may not be identified; check the specification.
#> [1] FALSE

Note, that instead of just assigning m to the second model m2 we need to clone() the object. R6 objects have reference semantics, so a direct assignment would refer to the same object. Cloning creates an independent object.