Marginal Effects of Mlogit in R

How to compute marginal effects of a multinomial logit model created with the nnet package?

The marginaleffects package should work in theory, but my example doesn't compile because of file size restrictions (meaning I don't have enough RAM for the 1.5 GB vector it tries to use). It's not even that large of a dataset, which is odd.

If you use marginal_effects() (margins package) for multinomial models, it only displays the output for a default category. You have to manually set each category you want to see. You can clean up the output with broom and then combine some other way. It's clunky, but it can work.

marginal_effects(model, category = 'cat1')

Marginal effects from the multinomial model

I believe you are missing the covariate information that needs to be put there, so if you use effects(mlfit, covariate = 'D'), It should work. Now the error is coming because the default of covariate is NULL. NULL is special in R, it has no(zero) length and hence you are getting argument of length zero. Please let me know if it fixes your issue.

As per documentation of effects.mlogit , it says:

covariate 
the name of the covariate for which the effect should be computed,

I am getting this output at my end:

R>effects(mlfit, covariate = 'D')
              1               2               3 
-0.003585105992 -0.070921137682 -0.026032167377 
              4               5 
 0.078295227196  0.022243183855

How to get marginal effects for categorical variables in mlogit?

It is kind of expected that effects doesn't work with factors since otherwise the output would contain another dimension, somewhat complicating the results, and it is quite reasonable that, just like in my solution below, one may instead want effects only for a certain factor level, rather than all the levels. Also, as I explain below, the marginal effects in the case of categorical variables are not uniquely defined, so that would be an additional complication for effects.

A natural workaround is to manually convert the factor variables into a series of dummy variables as in

aux <- model.matrix(y ~ x, df3F)[, -1]
head(aux)
#   x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13
# 1  0  0  0  0  0  0  0  0   0   1   0   0
# 2  0  0  0  0  0  0  0  0   0   1   0   0
# 3  0  0  0  0  0  1  0  0   0   0   0   0
# 4  0  0  0  0  0  0  0  0   1   0   0   0
# 5  0  0  0  0  0  0  0  1   0   0   0   0
# 6  0  0  0  0  0  0  1  0   0   0   0   0

so that the data then is

ml.datF <- mlogit.data(data.frame(y = df3F$y, aux), choice = "y", shape = "wide")

We also need to construct the formula manually with

frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse = " + "))))

So far so good. Now if we run

fit.mnlF <- mlogit(frml, data = ml.datF)
head(effects(fit.mnlF, covariate = "x2", data = ml.datF))
#           FALSE         TRUE
# 1 -1.618544e-15 0.000000e+00
# 2 -1.618544e-15 0.000000e+00
# 3 -7.220891e-08 7.221446e-08
# 4 -1.618544e-15 0.000000e+00
# 5 -5.881129e-08 5.880851e-08
# 6 -8.293366e-08 8.293366e-08

then the results are not correct. What effects did here is that it saw x2 as a continuous variable and computed the usual marginal effect for those cases. Namely, if the coefficient corresponding to x2 is b2 and our model is f(x,b2), effects computed the derivative of f with respect to b2 and evaluated at each observed vector x_i. This is wrong because x2 only takes values 0 and 1, not something around 0 or around 1, which is what taking the derivate assumes (the concept of a limit)! For instance, consider your other dataset df1. In that case we incorrectly get

colMeans(effects(fit.mnlF, covariate = "x2", data = ml.datF))
#           1           2           3           4           5 
# -0.25258378  0.07364406  0.05336283  0.07893391  0.04664298

Here's another way (using derivative approximation) to get this incorrect result:

temp <- ml.datF
temp$x2 <- temp$x2 + 0.0001
colMeans(predict(fit.mnlF, newdata = temp, type = "probabilities") - 
             predict(fit.mnlF, newdata = ml.datF, type = "probabilities")) / 0.0001
#           1           2           3           4           5 
# -0.25257597  0.07364089  0.05336032  0.07893273  0.04664202

Instead of using effects, I computed the wrong marginal effects by hand by using predict twice: the result is mean({fitted probability with x2new = x2old + 0.0001} - {fitted probability with x2new = x2old}) / 0.0001. That is, we looked at the change in the predicted probability by moving x2 up by 0.0001, which is either from 0 to 0.0001 or from 1 to 0.0001. Both of those don't make sense. Of course, we shouldn't expect anything else from effects since x2 in the data is numeric.

So then the question is how to compute right (average) marginal effects. As I said, marginal effect for categorical variables is not uniquely defined. Suppose that x_i is whether an individual i has a job and y_i is whether they have a car. So, then there are at least the following six things to consider.

Effect on the probability of y_i = 1 when going from x_i=0 to x_i=1.
When going from x_i=0 to x_i (the observed value).
From x_i to 1.

Now when we are interested in average marginal effects, we may want to average only over those individuals for whom the change in 1-3 makes a difference. That is,

From x_i=0 to x_i=1 if the observed value isn't 1.
From x_i=0 to x_i if the observed value isn't 0.
From x_i to 1 if the observed value isn't 1.

According to your results, Stata uses option 5, so I'll reproduce the same results, but it is straightforward to implement any other option, and I suggest to think which ones are interesting in your particular application.

AME.fun2 <- function(betas) {
  aux <- model.matrix(y ~ x, df1)[, -1]
  ml.datF <- mlogit.data(data.frame(y = df1$y, aux), choice="y", shape="wide")
  frml <- mFormula(formula(paste("y ~ 1 |", paste(colnames(aux), collapse=" + "))))
  fit.mnlF <- mlogit(frml, data = ml.datF)
  fit.mnlF$coefficients <- betas
  aux <- ml.datF # Auxiliary dataset
  aux$x3 <- 0 # Going from 0 to the observed x_i
  idx <- unique(aux[aux$x3 != ml.datF$x3, "chid"]) # Where does it make a change?
  actual <- predict(fit.mnlF, newdata = ml.datF)
  counterfactual <- predict(fit.mnlF, newdata = aux)
  colMeans(actual[idx, ] - counterfactual[idx, ])
}
(AME.mnl <- AME.fun2(ml.fit$coefficients))
#           1           2           3           4           5 
# -0.50000000 -0.17857142  0.03896104  0.16558441  0.47402597 

require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), nrow = 1, byrow = TRUE)
AME.mnl / AME.mnl.se
#           [,1]      [,2]    [,3]     [,4]     [,5]
# [1,] -5.291503 -2.467176 0.36922 1.485058 4.058994

Calculating marginal effects using mlogit in R

I sort of got an answer - not sure if its the best, but it worked. My covariate that I was interested in just so happened to be 2 classes - which means I could turn the covariate into a binary 0,1 numeric response. Therefore, when I rerun the code, I could then calculate the mean for this "categorical" variable.

However, I think that I am missing the point with this marginal effect and am likely using it or trying to interpret it incorrectly.

How to get average marginal effects (AMEs) with standard errors of a multinomial logit model?

We can do something very similar to what is done in your linked answer. In particular, first we want a function that would compute AMEs at a given vector of coefficients. For that we can define

AME.fun <- function(betas) {
  tmp <- ml.fit
  tmp$coefficients <- betas
  ME.mnl <- sapply(c.names, function(x) 
    effects(tmp, covariate = x, data = ml.d), simplify = FALSE)
  c(sapply(ME.mnl, colMeans))
}

where the second half is yours, while in the first one I use a trick to take the same ml.fit object and to change its coefficients. Next we find the jacobian with

require(numDeriv)
grad <- jacobian(AME.fun, ml.fit$coef)

and apply the delta method. Square roots of the diagonal of grad %*% vcov(ml.fit) %*% t(grad) is what we want. Hence,

(AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), nrow = 3, byrow = TRUE))
#             [,1]        [,2]        [,3]        [,4]        [,5]
# [1,] 0.003269320 0.004788536 0.004995723 0.004009762 0.002527462
# [2,] 0.004375795 0.006348496 0.008168883 0.008844684 0.005763966
# [3,] 0.009233616 0.014048212 0.014713090 0.012702188 0.008261734
AME.mnl / AME.mnl.se
#            1          2          3         4        5
# D  -9.259050 -1.8389907 0.30847523 4.2861720 8.051269
# x1 -6.657611 -2.4808393 1.59847852 1.4969683 3.224159
# x2 -2.950794 -0.3902812 0.05828811 0.4197057 3.212458

which coincides with Stata's results.