calculate marginal effects using mlogit package

lmoulin · Jun 7, 2013

Hi R-users

I try to calculate marginal effects of a multinomial logistic regression. To do this i use mlogit package and effects() function.

Here is how the procedure works (source : effects() function of mlogit package) :

Code:

data("Fishing", package = "mlogit")
Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
m <- mlogit(mode ~ price | income | catch, data = Fish)
# compute a data.frame containing the mean value of the covariates in the sample
z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean), 
catch = tapply(catch, index(m)$alt, mean), 
income = mean(income)))
# compute the marginal effects (the second one is an elasticity
effects(m, covariate = "income", data = z)
effects(m, covariate = "price", type = "rr", data = z)
effects(m, covariate = "catch", type = "ar", data = z)

I have no problem with first step (mlogit.data() function). I think my problem is on the specification of the multinomial regression.

My regression (for example with three variables) is on the form:

Code:

Y ~ 0 | X1 + X2 + X3

. When I try to estimate the marginal effects for model with 2 variables no problem, however for 3 variables R console returns me the following error: "Error in if (rhs% in% c (1, 3)) {: argument is of length zero " (translation from error in R console in french).

To understand what is my problem I tried to perform a multinomial regression of similar shape on the dataset "Fishing", that is to say :

Code:

mode ~ 0 | income + price + catch

(even if this form has no "economic" sense.) Again the R console returns me the same error for 3 variables but manages to estimate these effects for a model with two variables ...

This leads me to think that my problem really comes from the specification of my multinomial regression ... Do you know how I could find a solution to my problem?

Thank you for your help

chetan.apa · Jun 7, 2013

@lmoulin..Strangely, when I tried running your code and I didn't encounter any error. Below is the output I got

Code:

> data("Fishing", package = "mlogit")
> library(mlogit)
> Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
> m <- mlogit(mode ~ price | income | catch, data = Fish)
> # compute a data.frame containing the mean value of the covariates in the sample
> z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean), 
+                            catch = tapply(catch, index(m)$alt, mean), 
+                            income = mean(income)))
> # compute the marginal effects (the second one is an elasticity
> effects(m, covariate = "income", data = z)
        beach          boat       charter          pier 
 1.132964e-06  3.113068e-05 -2.408677e-05 -8.176877e-06 
> effects(m, covariate = "price", type = "rr", data = z)
             beach       boat    charter       pier
beach   -2.4634447  0.1512128  0.1512128  0.1512128
boat     0.5797049 -0.8172610  0.5797049  0.5797049
charter  0.9741351  0.9741351 -1.1590941  0.9741351
pier     0.1844507  0.1844507  0.1844507 -2.4302068
> effects(m, covariate = "catch", type = "ar", data = z)
               beach        boat     charter         pier
beach    0.040943135 -0.01803326 -0.01984425 -0.003065622
boat    -0.010447074  0.10568073 -0.08249023 -0.012743428
charter -0.012623645 -0.09057996  0.11860203 -0.015398428
pier    -0.001887074 -0.01354055 -0.01490036  0.030327981

lmoulin · Jun 7, 2013

Dear chetan.apa,

Than you for your help

but it's not exactly that i want.

I try to estimate the marginal effects for model of the form :

Code:

 Y ~ 0 | X1 + X2 + X3

Apply to Fishind data, m becomes

Code:

m <- mlogit(mode ~ 0 | price + income + catch, data = Fish)

.

But with this formula i can't have marginal effects and i don't understand why ...

chetan.apa · Jun 7, 2013

Ohh ok.. i looked at it again..I don't have much clue but it may have either something to do with the syntax of the model specified or probably that for this type of model marginal effects calculation would be different (i am just guessing here..)
someone with knowledge on this topic would be able to guide you better.
However I have an observation which might/might not help..the effects work for the last variable in the mlogit statement..weird observation I know

See code below

Code:

> data("Fishing", package = "mlogit")
> library(mlogit)
> Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
> #m <- mlogit(mode ~ price | income | catch, data = Fish)
> m <- mlogit(mode ~ 0 | catch + price + income, data = Fish)
> # compute a data.frame containing the mean value of the covariates in the sample
> z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean), 
+                            catch = tapply(catch, index(m)$alt, mean), 
+                            income = mean(income)))
> # compute the marginal effects (the second one is an elasticity
> effects(m, covariate = "income", data = z)
        beach          boat       charter          pier 
-1.942632e-05  4.328966e-05 -2.086361e-05 -2.999728e-06 
> effects(m, covariate = "price", type = "rr", data = z)
Error in if (rhs %in% c(1, 3)) { : argument is of length zero
> effects(m, covariate = "catch", type = "ar", data = z)
Error in if (rhs %in% c(1, 3)) { : argument is of length zero

Code:

> data("Fishing", package = "mlogit")
> library(mlogit)
> Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")
> #m <- mlogit(mode ~ price | income | catch, data = Fish)
> m <- mlogit(mode ~ 0 | catch + income +price, data = Fish)
> # compute a data.frame containing the mean value of the covariates in the sample
> z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean), 
+                            catch = tapply(catch, index(m)$alt, mean), 
+                            income = mean(income)))
> # compute the marginal effects (the second one is an elasticity
> effects(m, covariate = "income", data = z)
Error in if (rhs %in% c(1, 3)) { : argument is of length zero
> effects(m, covariate = "price", type = "rr", data = z)
     beach       boat    charter       pier 
 1.3373414 -0.2127258  0.1133452 -1.3069955 
> effects(m, covariate = "catch", type = "ar", data = z)
Error in if (rhs %in% c(1, 3)) { : argument is of length zero

lmoulin · Jun 7, 2013

I observed the same effect and I do not understand why only the latter effect each time is taken into account .. It's strange because the model I'm trying to estimate is the standard (see http://www.polsci.ucsb.edu/faculty/glasgow/ps206/ps206_mnl.r for example) multinomial logit model .. (with form Y ~ 0 | X1 + X2 + X3)

LauraS. · Aug 1, 2013

Hi Lmoulin,
I had the same problem with the effects comand a while ago. I solved it with a little transformation in the source code of effects.mlogit. In line 16 you should replace "cov.list <- lapply(attr(formula(object), "rhs"), as.character)" with "cov.list <- strsplit(as.character(attr(formula(object), "rhs")), " + ", fixed = TRUE)".
That should solve the problem.
best, laura

lmoulin · Aug 1, 2013

Hi Laura,

Thank you for your answer !!

Léonard.

Huckleberry · Jul 29, 2014

Hi Laura, I too am having the same problem on how to calculate marginal effects using mlogit package. But your wonderful answer has made me jump for joy - until I realized I have no idea how to do this "simple" solution. Can you provide details of how to do this? I use R studio, but I could easily switch if its necessary. Ive tried to use debug() but that gets me no where because effects.mlogit is actually using effects from stats library - and I just have no idea. ~help please

woudy12 · Aug 7, 2014

Hi Huckleberry,

I had the very same problem as you are describing here.. But after a daunting search i found the (not the most elegant i suppose) solution. I hope its not too late for you... In order to find the source code (it is apparently an S3 method) I used:

Code:

> methods(effects)
[1] effects.glm*    effects.lm*     effects.mlogit*

   Non-visible functions are asterisked
> capture.output(getAnywhere(effects.mlogit), file='Your Desired Filename.r')

The source code of effects.mlogit is now captured in the r file. When you open it, give a (random) object name to the function by using "Obj Name <- Function.... " in the first line of the code, replace the 16th line as Laura stated (with "cov.list...") and save the file it should work. The function is then called by the Obj name you've given to it, after loading your r file by the following code:

Code:

> source(file="Your Desired Filename.r")
>#effects(mlogit."Your Model", covariate = "Your Covariate", data = "Your Data")  becomes
>"Obj Name"(mlogit."Your Model", covariate = "Your Covariate", data = "Your Data")

I am using Rstudio btw..
Hope this helps for you 2..

Cheers, Wouter

NB: Thanks to this post http://stackoverflow.com/questions/19226816/how-can-i-view-the-source-code-for-a-function

Thomas.V · Mar 12, 2016

For those who are still unable, such as me, to follow the steps suggested by Laura or Wouter, simply assign the repaired effects.mlogit function - thanks to Laura - to a new name, say m.effects. This you do by running the code below. Then you can simply use m.effects(...) instead of effects(...).

Code:

m.effects <- function (object, covariate = NULL, type = c("aa", "ar", "rr", 
    "ra"), data = NULL, ...) 
{
    type <- match.arg(type)
    if (is.null(data)) {
        P <- predict(object, returnData = TRUE)
        data <- attr(P, "data")
        attr(P, "data") <- NULL
    }
    else P <- predict(object, data)
    newdata <- data
    J <- length(P)
    alt.levels <- names(P)
    pVar <- substr(type, 1, 1)
    xVar <- substr(type, 2, 2)
    cov.list <- strsplit(as.character(attr(formula(object), "rhs")), " + ", fixed = TRUE)
    rhs <- sapply(cov.list, function(x) length(na.omit(match(x, 
        covariate))) > 0)
    rhs <- (1:length(cov.list))[rhs]
    eps <- 1e-05
    if (rhs %in% c(1, 3)) {
        if (rhs == 3) {
            theCoef <- paste(alt.levels, covariate, sep = ":")
            theCoef <- coef(object)[theCoef]
        }
        else theCoef <- coef(object)[covariate]
        me <- c()
        for (l in 1:J) {
            newdata[l, covariate] <- data[l, covariate] + eps
            newP <- predict(object, newdata)
            me <- rbind(me, (newP - P)/eps)
            newdata <- data
        }
        if (pVar == "r") 
            me <- t(t(me)/P)
        if (xVar == "r") 
            me <- me * matrix(rep(data[[covariate]], J), J)
        dimnames(me) <- list(alt.levels, alt.levels)
    }
    if (rhs == 2) {
        newdata[, covariate] <- data[, covariate] + eps
        newP <- predict(object, newdata)
        me <- (newP - P)/eps
        if (pVar == "r") 
            me <- me/P
        if (xVar == "r") 
            me <- me * data[[covariate]]
        names(me) <- alt.levels
    }
    me
}

Hope that helps.
Thomas

calculate marginal effects using mlogit package

lmoulin

New Member

chetan.apa

Member

lmoulin

New Member

chetan.apa

Member

lmoulin

New Member

LauraS.

New Member

lmoulin

New Member

Huckleberry

New Member

woudy12

New Member

Thomas.V

New Member