Странное поведение команды эффектов mlogit в R

Я оцениваю мультиномиальную логит-модель и хотел бы сообщить о предельном эффекте. Я столкнулся с трудностью, так как при использовании более крупной версии модели я получаю сообщение об ошибке.

Вот воспроизводимый пример. Следующий код с двумя ковариатами работает нормально.

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)

Теперь, когда у меня есть три ковариаты:

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)

Последняя строка выдает следующую ошибку:

Ошибка в if (rhs %in% c(1, 3)) { : аргумент имеет нулевую длину

Но если я побегу

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

тогда он работает нормально для получения предельных эффектов col3. Почему бы ему не дать предельных эффектов col1?

Обратите внимание, что все столбцы не содержат NULL и имеют одинаковую длину. Может кто-нибудь объяснить, в чем причина такого поведения?


r regression logistic-regression mlogit non-linear-regression
person splinter    schedule 11.04.2017    source источник
comment
Это означает, что rhs — пустой вектор. Вам, вероятно, придется сделать это воспроизводимым, чтобы получить полный ответ. Один маленький шаг к воспроизводимости — это отметить, какая строка вызывает ошибку.   -  person lmo    schedule 11.04.2017
comment
Это последняя строка, я упоминал об этом в вопросе (в конце)   -  person splinter    schedule 11.04.2017
comment
Одна из возможностей — покопаться в функции effects (если она не скомпилирована), чтобы увидеть, на что ссылается rhs. Тогда вы сможете вернуться к тому, почему он пуст.   -  person lmo    schedule 11.04.2017
comment
воспроизводимый пример значительно повысит ваши шансы на получение ответа (кто-то другой может захотеть покопаться в effects за вас, но если вы не приведете воспроизводимый пример, им придется придумать свой собственный воспроизводимый пример, который будет намного сложнее и, возможно, не стоит)   -  person Ben Bolker    schedule 11.04.2017
comment
Спасибо @BenBolker! Добавлен воспроизводимый пример.   -  person splinter    schedule 11.04.2017
comment
@splinter - указатель и код решили вашу проблему? Я считаю, что это так - просто проверяю, не пропустил ли я что-нибудь.   -  person Technophobe01    schedule 16.04.2017
comment
Большое спасибо, @Technophobe01, я пока не могу проверить это, потому что у меня нет другого компьютера, но я получу сдачу в ближайшие несколько дней.   -  person splinter    schedule 16.04.2017
comment
@splinter - Не беспокойтесь - рад помочь. Шотландия место находки. Посмотрите на Глена Коу - я скучаю по пейзажам :-)   -  person Technophobe01    schedule 16.04.2017
comment
Спасибо @Technophobe01, Глен Коу очень мил :)   -  person splinter    schedule 19.04.2017

Ответы (1)


arrow_upward
1
arrow_downward

Я чувствую, что это может помочь вам найти решение.

Ссылка: 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

Объяснение:

Требуется небольшое преобразование в исходном коде Effects.mlogit.

В строке 16 вы должны заменить "cov.list ‹- lapply(attr(formula(object), "rhs"), as.character)" на "cov.list ‹- strsplit(as.character(attr(formula(object)) , "правая")), "+", фиксированное = ИСТИНА)"

Исправить результат:

> 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

Код

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)
person Technophobe01    schedule 15.04.2017