Strange behavior of the effects command of mlogit in R

Question

I am estimating a multionomial logit model and would like to report marginal effect. I am running into a difficulty, as when I am using a larger version of the model I get an error.

Here is an reproducible example. The following code, with two covariates, works fine.

library(mlogit)

df = data.frame(c(0,1,1,2,0,1,0), c(1,6,7,4,2,2,1), c(683,276,756,487,776,100,982))
colnames(df) <- c('y', 'col1', 'col3')
df$col2<-df$col1^2
mydata = df

mldata <- mlogit.data(mydata, choice="y", shape="wide")
mlogit.model1 <- mlogit(y ~ 1| col1+col2, data=mldata)
m <- mlogit(y ~ 1| col1+col2, data = mldata)
z <- with(mldata, data.frame(col1 = tapply(col1, index(m)$alt, mean), 
                             col2 = tapply(col2, index(m)$alt, mean) ) )
effects(mlogit.model1, covariate = "col1", data = z)

Now, when I have three covariates:

mlogit.model1 <- mlogit(y ~ 1| col1+col2+col3, data=mldata)
m <- mlogit(y ~ 1| col1+col2+col3, data = mldata)
z <- with(mldata, data.frame(col1 = tapply(col1, index(m)$alt, mean), 
                             col2 = tapply(col2, index(m)$alt, mean), 
                             col3 = tapply(col3, index(m)$alt, mean) ) )
effects(mlogit.model1, covariate = "col1", data = z)

The last line gives the following error:

Error in if (rhs %in% c(1, 3)) { : argument is of length zero

But if I run

effects(mlogit.model1, covariate = "col3", data = z)

then it works ok for giving the marginal effects of col3. Why would it not give the marginal effects of col1?

Note that all columns contain no NULLs and are of the same length. Can someone explain what's the reason for this behavior?

Answer 1

My sense is this may help guide you to a solution.

Reference: http://www.talkstats.com/showthread.php/44314-calculate-marginal-effects-using-mlogit-package

> methods(effects)
[1] effects.glm*    effects.lm*     effects.mlogit*
see '?methods' for accessing help and source code 
Note: Non-visible functions are asterisked

Explanation:

A little transformation in the source code of effects.mlogit is required.

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)"

Fix Result:

> effects(mlogit.model1, covariate = "col1", data = z)
            0             1             2 
-4.135459e-01  4.135459e-01  9.958986e-12 

> myeffects(mlogit.model2, covariate = "col1", data = z2)
           0            1            2 
 1.156729129 -1.157014778  0.000285649 

Code

require(mlogit)

myeffects<-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
}

df = data.frame(c(0,1,1,2,0,1,0), c(1,6,7,4,2,2,1), c(683,276,756,487,776,100,982))
colnames(df) <- c('y', 'col1', 'col3')
df$col2<-df$col1^2
mydata = df

mldata <- mlogit.data(mydata, choice="y", shape="wide")
mlogit.model1 <- mlogit(y ~ 1| col1+col2, data=mldata)
m <- mlogit(y ~ 1| col1+col2, data = mldata)
z <- with(mldata, data.frame(col1 = tapply(col1, index(m)$alt, mean), 
                             col2 = tapply(col2, index(m)$alt, mean) ) )

mldata2 <- mlogit.data(mydata, choice="y", shape="wide")
mlogit.model2 <- mlogit(y ~ 1| col1+col2+col3, data=mldata2)
m2 <- mlogit(y ~ 1| col1+col2+col3, data = mldata2)
z2 <- with(mldata, data.frame(col1 = tapply(col1, index(m2)$alt, mean), 
                             col2 = tapply(col2, index(m2)$alt, mean), 
                             col3 = tapply(col3, index(m2)$alt, mean) ) )

effects(mlogit.model1, covariate = "col1", data = z)
myeffects(mlogit.model2, covariate = "col1", data = z2)