Where is the variance-covariance matrix for the coefficients from the command mlogit from the R package mlogit stored or how can I calculate it?

Question

I have estimated a nested logit using the command mlogit() from the R package mlogit. The output is in the form of a list with 17 elements. If I run summary() on the output, I get a table with the estimated coefficients and their respective standard errors. However, I cannot find such standard errors in the output of the mlogit() command (the list with 17 elements). Am I missing where it is? Below I provide a non-reproducible example of the kind of code I am running, as well as the list of elements in the output.

My guess is that maybe the command summary() is somehow calculating the standard errors for the estimated coefficients. Given the coefficients are estimated by mlogit() using Maximum Likelihood, I am guessing these standard errors come from calculating the inverse of the Fischer information matrix. However, I am unsure how to do so given the output I have from the mlogit() command. I will be thankful if anyone could help me with this.

Thanks in advance for your help!

CODE EXAMPLE:

library(dfidx)
library(mlogit)

ml_data <- data %>% dfidx(idx = c("id", "alternative"))
model_output <- mlogit(choice ~ regresor1 + regresor2,
                   data = ml_data,
                   reflevel = "alt1",
                   nests = list(nest1 = c("alt2", "alt3", "alt4"),
                                nest2 = c("alt1")),
                   un.nest.el = FALSE, constPar = c("iv:nest2" = 1))
summary(model_output)

OUTPUT FROM mlogit():

More information available here: Rdocumentation (though my list has 17 elements and the documentation only describes 13 of them)

See here: https://rdrr.io/cran/mlogit/src/R/methods.R for the details of `summary.mlogit()` — Paul Stafford Allen, Jul 07 '23 at 13:27

score 2 · Answer 1 · answered Jul 07 '23 at 13:27

see ?vcov.mlogit

you can also inspect the function by typing mlogit:::summary.mlogit and hitting ENTER (note the triple colon to expose the "hidden" machinery of an R package)

You will find both options by browsing the suggestions popping up after hitting TAB twice. If you already have an idea of what you're looking for (e. g. summary.xyz or plot.xyz for a method specific to a result of class xyz), type the known fragment to narrow down the results of TAB expansion.

score 1 · Accepted Answer · answered Jul 07 '23 at 13:36

Your example isn't reproducible, but we can run the example from ?mlogit:

library(mlogit)

data("Fishing", package = "mlogit")
Fish <- dfidx(Fishing, varying = 2:9, shape = "wide", choice = "mode")

model <- mlogit(mode ~ price + catch, data = Fish)

This results in

model
#> 
#> Call:
#> mlogit(formula = mode ~ price + catch, data = Fish, method = "nr")
#> 
#> Coefficients:
#>    (Intercept):boat  (Intercept):charter     (Intercept):pier  
#>             0.87137              1.49889              0.30706  
#>               price                catch  
#>            -0.02479              0.37717

Notice that the class of the model object is "mlogit":

class(model)
#> [1] "mlogit"

In R, the summary function is an S3 generic, and when we call summary(model), it actually calls the mlogit:::summary.mlogit function:

summary(model)
#> 
#> Call:
#> mlogit(formula = mode ~ price + catch, data = Fish, method = "nr")
#> 
#> Frequencies of alternatives:choice
#>   beach    boat charter    pier 
#> 0.11337 0.35364 0.38240 0.15059 
#> 
#> nr method
#> 7 iterations, 0h:0m:0s 
#> g'(-H)^-1g = 6.22E-06 
#> successive function values within tolerance limits 
#> 
#> Coefficients :
#>                       Estimate Std. Error  z-value  Pr(>|z|)    
#> (Intercept):boat     0.8713749  0.1140428   7.6408 2.154e-14 ***
#> (Intercept):charter  1.4988884  0.1329328  11.2755 < 2.2e-16 ***
#> (Intercept):pier     0.3070552  0.1145738   2.6800 0.0073627 ** 
#> price               -0.0247896  0.0017044 -14.5444 < 2.2e-16 ***
#> catch                0.3771689  0.1099707   3.4297 0.0006042 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Log-Likelihood: -1230.8
#> McFadden R^2:  0.17823 
#> Likelihood ratio test : chisq = 533.88 (p.value = < 2.22e-16)

This doesn't simply print a summary to the screen - it actually does further calculations that produce an object of class "summary.mlogit". We can see this by storing the object and examining it:

model_summary <- summary(model)

class(model_summary)
#> [1] "summary.mlogit" "mlogit"

We will see that it contains different data members than model does:

names(model_summary)
#>  [1] "coefficients"  "logLik"        "gradient"      "hessian"      
#>  [5] "est.stat"      "fitted.values" "probabilities" "linpred"      
#>  [9] "indpar"        "residuals"     "omega"         "rpar"         
#> [13] "nests"         "model"         "freq"          "formula"      
#> [17] "call"          "CoefTable"     "lratio"        "mfR2"

One of these, the member CoefTable is the data table printed in the summary:

model_summary$CoefTable
#>                        Estimate  Std. Error    z-value     Pr(>|z|)
#> (Intercept):boat     0.87137491 0.114042831   7.640769 2.153833e-14
#> (Intercept):charter  1.49888838 0.132932796  11.275535 0.000000e+00
#> (Intercept):pier     0.30705525 0.114573796   2.679978 7.362701e-03
#> price               -0.02478955 0.001704403 -14.544420 0.000000e+00
#> catch                0.37716885 0.109970659   3.429723 6.041986e-04

If we don't want to look through all the members of an object to find the summary table, we can always just try the generic S3 coef, to see if the authors have defined a method. It turns out they have:

coef(model_summary)
#>                        Estimate  Std. Error    z-value     Pr(>|z|)
#> (Intercept):boat     0.87137491 0.114042831   7.640769 2.153833e-14
#> (Intercept):charter  1.49888838 0.132932796  11.275535 0.000000e+00
#> (Intercept):pier     0.30705525 0.114573796   2.679978 7.362701e-03
#> price               -0.02478955 0.001704403 -14.544420 0.000000e+00
#> catch                0.37716885 0.109970659   3.429723 6.041986e-04

It seems that this is what you are looking for. However, you mention a variance-covariance matrix in your question, which is a different concept from the coefficient table. If you genuinely want the variance-covariance matrix, there is a vcov method that you can use on your model too:

vcov(model)
#>                     (Intercept):boat (Intercept):charter (Intercept):pier
#> (Intercept):boat        1.300577e-02        1.025681e-02     7.519846e-03
#> (Intercept):charter     1.025681e-02        1.767113e-02     7.145189e-03
#> (Intercept):pier        7.519846e-03        7.145189e-03     1.312715e-02
#> price                  -1.807016e-06       -9.194307e-05     8.030833e-07
#> catch                   7.717477e-04       -5.247810e-03     7.401401e-04
#>                             price         catch
#> (Intercept):boat    -1.807016e-06  7.717477e-04
#> (Intercept):charter -9.194307e-05 -5.247810e-03
#> (Intercept):pier     8.030833e-07  7.401401e-04
#> price                2.904989e-06  1.010511e-05
#> catch                1.010511e-05  1.209355e-02

^{Created on 2023-07-07 with reprex v2.0.2}

Where is the variance-covariance matrix for the coefficients from the command mlogit from the R package mlogit stored or how can I calculate it?

2 Answers2