Estimate multinomial probit model with mlogit (R package)

From the document and help, probit model is supported by mlogit. But when I tried it with these R scripts, the estimation takes much longer time to run (than the logit verion) and the result is quite a bit different as well (argument probit=FALSE). Does the probit behave correctly? If so, how should I interpret the er.gc, er.gr, etc coefficients?

> require(mlogit)
> data(Heating)
> H <- mlogit.data(Heating, shape="wide", choice="depvar", varying=c(3:12))
> m1.probit = mlogit(depvar~ic+oc, H, probit=TRUE)
> summary(m1.probit)

Call:
mlogit(formula = depvar ~ ic + oc, data = H, probit = TRUE)

Frequencies of alternatives:
ec       er       gc       gr       hp
0.071111 0.093333 0.636667 0.143333 0.055556

bfgs method
37 iterations, 0h:4m:54s
g'(-H)^-1g = 0.011
last step couldn't find higher value

Coefficients :
Estimate  Std. Error t-value Pr(>|t|)
er:(intercept)  2.5611e-01  3.6641e-01  0.6990  0.48457
gc:(intercept) -2.6944e-02  3.3211e-01 -0.0811  0.93534
gr:(intercept) -1.8439e+01  3.2798e+01 -0.5622  0.57398
hp:(intercept) -6.4231e-01  7.4214e-01 -0.8655  0.38677
ic             -1.1447e-03  5.3175e-04 -2.1528  0.03133 *
oc             -3.3779e-03  1.4011e-03 -2.4109  0.01591 *
er.gc           4.4987e-01  2.6880e-01  1.6736  0.09421 .
er.gr           5.8580e+00  1.1236e+01  0.5214  0.60212
er.hp           1.2613e+00  5.0231e-01  2.5109  0.01204 *
gc.gc           7.1013e-01  3.5489e-01  2.0010  0.04540 *
gc.gr          -8.4606e+00  1.6848e+01 -0.5022  0.61555
gc.hp           6.7245e-01  6.1475e-01  1.0939  0.27401
gr.gr           1.4085e+01  2.6034e+01  0.5410  0.58849
gr.hp           4.9476e-01  4.4568e-01  1.1101  0.26694
hp.hp           2.2620e-01  2.9062e-01  0.7783  0.43637
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1000.1
Likelihood ratio test : chisq = 44.213 (p.value = 1.4049e-05)

> m1.logit = mlogit(depvar~ic+oc, H, probit=FALSE)
> summary(m1.logit)

Call:
mlogit(formula = depvar ~ ic + oc, data = H, probit = FALSE,
method = "nr", print.level = 0)

Frequencies of alternatives:
ec       er       gc       gr       hp
0.071111 0.093333 0.636667 0.143333 0.055556

nr method
6 iterations, 0h:0m:0s
g'(-H)^-1g = 9.58E-06
successive function values within tolerance limits

Coefficients :
Estimate  Std. Error t-value  Pr(>|t|)
er:(intercept)  0.19459102  0.20424212  0.9527 0.3407184
gc:(intercept)  0.05213336  0.46598878  0.1119 0.9109210
gr:(intercept) -1.35058266  0.50715442 -2.6631 0.0077434 **
hp:(intercept) -1.65884594  0.44841936 -3.6993 0.0002162 ***
ic             -0.00153315  0.00062086 -2.4694 0.0135333 *
oc             -0.00699637  0.00155408 -4.5019 6.734e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -1008.2
Likelihood ratio test : chisq = 27.99 (p.value = 8.3572e-07)

• You have a nice fast machine -- took my computer 2x as long to fit that probit model! – atiretoo Jun 3 '13 at 21:18

The run with probit=TRUE has not converged to a good answer. See the line in the output that starts with 'last step could not find higher value' and compare the same section in the logit model output. The other reason it takes so long to fit the probit model is that the software is approximating a high dimensional integral using simulation (See the vignette for mlogit, pg 54). Sometimes rescaling covariates can help with numerical difficulties but that's not the case here, I tried

HS$ic <- scale(H$ic)
HS$oc <- scale(H$oc)

m2.probit = mlogit(depvar~ic+oc, HS, probit=TRUE)


and had the same difficulty. I would treat the results of the probit model with a degree of skepticism. In particular, there is something funny going on with the outcome 'gr' in the probit model (see intercepts, and variance parameter estimates).

The coefficients in the summary that are labeled er.gc, er.gr etc. are the parameters of the variance-covariance matrix that is being estimated as part of the probit model.

Probit models often take longer to fit because the likelihood function is calculated by simulation or quadrature. The logit likelihood has a closed form solution that makes it fast.

Also, the Probit likelihood function is not globally convex, so the algorithm can converge to local maxima. You need to try different starting values.

Finally, the coefficients should not be the same, because they use a different scale parameter.

In general, if there is not a really compelling reason to use Probit, it's best to just stay away.