Как получить предельные эффекты для категориальных переменных в mlogit? - PullRequest
Новая
       98

Как получить предельные эффекты для категориальных переменных в mlogit?

0 голосов
/ 09 января 2019

Я хочу вычислить предельные эффекты для объекта "mlogit", в котором объясняющие переменные являются категориальными (факторами). В то время как с числовыми данными effects() что-то выбрасывает, с категориальными данными это не так.

Для простоты ниже приведен двумерный пример.

числовые переменные 

# with mlogit
library(mlogit)
ml.dat <- mlogit.data(df3, choice="y", shape="wide")
fit.mnl <- mlogit(y ~ 1 | x, data=ml.dat)

head(effects(fit.mnl, covariate="x", data=ml.dat))
#         FALSE       TRUE
# 1 -0.01534581 0.01534581
# 2 -0.01534581 0.01534581
# 3 -0.20629452 0.20629452
# 4 -0.06903946 0.06903946
# 5 -0.24174312 0.24174312
# 6 -0.39306240 0.39306240

# with glm
fit.glm <- glm(y ~ x, df3, family = binomial)

head(effects(fit.glm))
# (Intercept)           x                                                 
#  -0.2992979  -4.8449254   2.3394989   0.2020127   0.4616640   1.0499595

факторные переменные 

# transform to factor
df3F <- within(df3, x <- factor(x))
class(df3F$x) == "factor"
# [1] TRUE

Пока glm() еще что-то выбрасывает, 

# with glm
fit.glmF <- glm(y ~ x, df3F, family = binomial)

head(effects(fit.glmF))
# (Intercept)           x2           x3           x4           x5           x6 
# 0.115076511 -0.002568206 -0.002568206 -0.003145397 -0.003631992 -0.006290794

подход mlogit() 

# with mlogit
ml.datF <- mlogit.data(df3F, choice="y", shape="wide")
fit.mnlF <- mlogit(y ~ 1 | x, data=ml.datF)

head(effects(fit.mnlF, covariate="x", data=ml.datF))

выдает это ошибка : 

Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels
In addition: Warning message:
In Ops.factor(data[, covariate], eps) :

 Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

Как я мог решить эту проблему?

Я уже пытался манипулировать effects.mlogit() с помощью этим ответом , но это не помогло решить мою проблему.

Примечание: Этот вопрос относится к этому решению , которое я хочу применить к категориальным объясняющим переменным.

редактировать

(Чтобы продемонстрировать проблему при применении данного решения к основной проблеме, связанной с вопросом, связанным выше. См. Комментарии.) 

# new example ----
library(mlogit)
ml.d <- mlogit.data(df1, choice="y", shape="wide")
ml.fit <- mlogit(y ~ 1 | factor(x), reflevel="1", data=ml.d)

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  # probably?
  colMeans(effects(fit.mnlF, covariate="x2", data=ml.datF))  # first co-factor?
}

(AME.mnl <- AME.fun2(ml.fit$coefficients))

require(numDeriv)
grad <- jacobian(AME.fun2, ml.fit$coef)
(AME.mnl.se <- matrix(sqrt(diag(grad %*% vcov(ml.fit) %*% t(grad))), 
                      nrow=3, byrow=TRUE))
AME.mnl / AME.mnl.se
#  doesn't work yet though...

# probably "true" values, obtained from Stata:
# # ame
#         1      2      3      4      5
# 1.     NA     NA     NA     NA     NA   
# 2. -0.400  0.121 0.0971  0.113 0.0686   
# 3. -0.500 -0.179 0.0390  0.166 0.474 
#
# # z-values
#        1     2     3     4     5
# 1.    NA    NA    NA    NA    NA
# 2. -3.86  1.25  1.08  1.36  0.99
# 3. -5.29 -2.47  0.37  1.49  4.06

данные 

df3 <- structure(list(x = c(11, 11, 7, 10, 9, 8, 9, 6, 9, 9, 8, 9, 11, 
7, 8, 11, 12, 5, 8, 8, 11, 6, 13, 12, 5, 8, 7, 11, 8, 10, 9, 
10, 7, 9, 2, 10, 3, 6, 11, 9, 7, 8, 4, 12, 8, 12, 11, 9, 12, 
9, 7, 7, 7, 10, 4, 10, 9, 6, 7, 8, 9, 13, 10, 8, 10, 6, 7, 10, 
9, 6, 4, 6, 6, 8, 6, 9, 3, 7, 8, 2, 8, 6, 7, 9, 10, 8, 6, 5, 
5, 7, 9, 1, 6, 11, 11, 9, 7, 8, 9, 9), y = c(TRUE, TRUE, TRUE, 
TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, 
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, 
TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, 
TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, 
TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, 
TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, 
FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, 
FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, 
FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, 
TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE
)), class = "data.frame", row.names = c(NA, -100L))

> summary(df3)
       x             y          
 Min.   : 1.00   Mode :logical  
 1st Qu.: 7.00   FALSE:48       
 Median : 8.00   TRUE :52       
 Mean   : 8.08                  
 3rd Qu.:10.00                  
 Max.   :13.00  

df1 <- structure(list(y = c(5, 4, 2, 2, 2, 3, 5, 4, 1, 1, 2, 4, 1, 4, 
5, 5, 2, 3, 3, 5, 5, 3, 2, 4, 5, 1, 3, 3, 4, 3, 5, 2, 4, 4, 5, 
5, 5, 2, 1, 5, 1, 3, 1, 4, 1, 2, 2, 4, 3, 1, 4, 3, 1, 1, 5, 2, 
5, 4, 2, 2, 4, 2, 3, 5, 4, 1, 2, 2, 3, 5, 2, 5, 3, 3, 3, 1, 3, 
1, 1, 4, 3, 4, 5, 2, 1, 1, 3, 1, 5, 4, 4, 2, 5, 3, 4, 4, 3, 1, 
5, 2), x = structure(c(2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 
2L, 1L, 1L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 3L, 3L, 2L, 
3L, 2L, 2L, 2L, 3L, 2L, 1L, 3L, 2L, 3L, 3L, 1L, 1L, 3L, 2L, 2L, 
1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 3L, 2L, 
2L, 2L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 1L, 3L, 2L, 2L, 1L, 2L, 
2L, 1L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 1L, 1L, 1L, 3L, 2L, 
3L, 2L, 3L, 1L, 2L, 3L, 3L, 1L, 2L, 2L), .Label = c("1", "2", 
"3"), class = "factor")), row.names = c(NA, -100L), class = "data.frame")

1 Ответ

0 голосов
/ 09 января 2019

Отчасти ожидается, что effects не будет работать с факторами, поскольку в противном случае выходные данные будут содержать другое измерение, несколько усложняющее результаты, и вполне разумно, что, как и в моем решении ниже, можно вместо этого хотеть эффекты только для определенного факторного уровня, а не для всех уровней. Кроме того, как я объясню ниже, предельные эффекты в случае категориальных переменных не определены однозначно, поэтому это будет дополнительным осложнением для effects.

Естественным обходным путем является ручное преобразование факторных переменных в ряд фиктивных переменных, как в 

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

так, чтобы данные тогда были 

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

Нам также нужно построить формулу вручную с помощью 

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

Пока все хорошо. Теперь, если мы запустим 

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

тогда результаты не верны. effects сделал здесь то, что он увидел x2 как непрерывную переменную и вычислил обычный предельный эффект для этих случаев. А именно, если коэффициент, соответствующий x2, равен b2, а наша модель - f (x, b2), effects вычислил производную от f по b2 и вычислил для каждого наблюдаемого вектора x i . Это неверно, потому что x2 принимает только значения 0 и 1, а не что-то около 0 или около 1, что предполагает принятие производного (концепция предела)! Например, рассмотрим ваш другой набор данных df1. В этом случае мы неправильно получаем 

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

Вот еще один способ (используя производную аппроксимацию), чтобы получить этот неверный результат: 

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

Вместо использования effects я вручную вычислял неправильные предельные эффекты, дважды используя predict: результат является средним ({соответствует вероятность с x2new = x2old + 0.0001} - {соответствует вероятность с x2new = x2old}) / 0,0001. То есть мы смотрели на изменение прогнозируемой вероятности, перемещая x2 на 0,0001, что составляет либо от 0 до 0,0001, либо от 1 до 0,0001. Оба из них не имеют смысла. Конечно, от effects не следует ожидать ничего другого, поскольку x2 в данных является числовым.

Итак, вопрос в том, как рассчитать правильные (средние) предельные эффекты. Как я уже сказал, предельный эффект для категориальных переменных не определяется однозначно. Предположим, что x_i - это ли человек, у которого есть работа, а y_i, есть ли у него машина. Итак, есть как минимум следующие шесть вещей, которые следует учитывать.

  1. Влияние на вероятность y_i = 1 при переходе от x_i = 0 к x_i = 1.
  2. При переходе от x_i = 0 к x_i (наблюдаемое значение).
  3. От x_i до 1.

Теперь, когда мы заинтересованы в средних предельных эффектах, мы можем захотеть усреднить только по тем людям, для которых изменение в 1-3 имеет значение. То есть

  1. От x_i = 0 до x_i = 1, если наблюдаемое значение не равно 1.
  2. От x_i = 0 до x_i, если наблюдаемое значение не равно 0.
  3. От x_i до 1, если наблюдаемое значение не равно 1.

Согласно вашим результатам, Stata использует вариант 5, поэтому я воспроизведу те же результаты, но реализовать любой другой вариант несложно, и я предлагаю подумать, какие из них интересны для вашего конкретного приложения. 

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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...