How to use multinomial probit coefficients to predict?

Question

I fitted a multinomial probit model with one independent categorical variable Y (levels 1,2,3) and two explanatory variables X1 and X2.

Using mlogit package in R like this:

library(mlogit)
df = read.csv("https://gitlab.com/cristiandavidarteaga/rtraining/raw/b40daf27a52bf01ce58d0ea32c5e4854f5b23836/mlogit_2var/data.csv",header = T)
d = mlogit.data(df,shape = "wide",choice = "y")

myprobit = mlogit(y~0|x1+x2, d, probit = TRUE)
summary(myprobit)

Gives me the following coefficients:

Frequencies of alternatives:
    1     2     3 
0.509 0.128 0.363 

bfgs method
21 iterations, 0h:0m:34s 
g'(-H)^-1g = 9.56E-08 
gradient close to zero 

Coefficients :
                 Estimate  Std. Error  t-value Pr(>|t|)    
2:(intercept) -10.7685665   0.9330425 -11.5413   <2e-16 ***
3:(intercept) -11.4357413   1.0913296 -10.4787   <2e-16 ***
2:x1            0.1097622   0.0093004  11.8019   <2e-16 ***
3:x1            0.1094478   0.0094566  11.5737   <2e-16 ***
2:x2            0.1010603   0.0100107  10.0952   <2e-16 ***
3:x2            0.1150660   0.0116610   9.8676   <2e-16 ***
2.3             0.9781048   0.0471720  20.7348   <2e-16 ***
3.3             0.0676135   0.0521005   1.2978   0.1944    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -199.84
McFadden R^2:  0.79498 
Likelihood ratio test : chisq = 1549.8 (p.value = < 2.22e-16)

I can't find a clear explanation about how to use these coefficients to predict outcomes for new data.

For example, If I have these coefficients, How can I manually predict (by hand, not using R) the outcome (1, 2 or 3) for x1 = 26 and x2 = 55 ?

Do I need to use the co-variance matrix to do this?

I know R or STATA can do it, however, for my research it's important to understand how to do it since I need to write a custom version of probit.

Answer 1

For example, If I have these coefficients, How can I manually predict (by hand, not using R) the outcome (1, 2 or 3) for x1 = 26 and x2 = 55 ?

You cannot do it by hand. The conditional probabilities are given by 2D integrals with the integrand being a multivariate normal distribution density function. There is no closed form solution. See the vignette called "6. The multinomial probit model" in version 1.1-1 of the mlogit package. You can also see my formulas in this question.

Computation of Multivariate Normal and t Probabilities by Genz and Bretz has a section on the 2D case but there is no closed form solution.

Do I need to use the co-variance matrix to do this?

Yes. See the aforementioned vignette.

Answer 2

In the Multinomial Probit model, recall that the probability of individual $i$ choosing alternative $j$, expressed as $\pi_{ij}$, is given by: $$P(Y = y_{ij}) = \pi_{ij} = \Phi(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \theta_j)$$ where $\Phi(z)$ = the standard normal cumulative density function, expressed as: $$\Phi(z) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{z}e^{-x^2/2}dx$$ Approximate values for $\Phi(z)$ can be found in tables in most statistics textbooks.

Hence from your results, the corresponding $\pi_{ij}$ are expressed as:

• For $j$ = 2 $$\pi_{i2} = \Phi(-10.769 + 0.11x_{1i} + 0.101x_{2i} + 0.978)$$
• For $j$ = 3 $$\pi_{i3} = \Phi(-11.436 + 0.109x_{1i} + 0.115x_{2i} + 0.068)$$