polr breaks with `unit_ses=TRUE`

[vincentarelbundock]

• [x] a possible bug
• [ ] a question about package functionality
• [ ] a suggested code or documentation change, improvement to the code, or feature request

library(MASS)
library(margins)

# Simulate dataset
n <- 100
u <- rnorm(n)
X1 <- u + rnorm(n)
X2 <- u + rnorm(n)
Ystar <- X1 + X2 + rnorm(n)
Y <- cut(Ystar, 3)
dat <- data.frame(X1, X2, Y, Ystar)

# Estimate
mod <- polr(Y ~ X1 + X2, data = dat)
margins(mod, unit_ses = FALSE)
#> 
#> Re-fitting to get Hessian
#> Average marginal effects
#> polr(formula = Y ~ X1 + X2, data = dat)
#>       X1      X2
#>  -0.1343 -0.1299
margins(mod, unit_ses = TRUE)
#> 
#> Re-fitting to get Hessian
#> Error in jacobian %*% vcov: non-conformable arguments

^{Created on 2019-07-16 by the reprex package (v0.3.0)}

> traceback()
7: diag(jacobian %*% vcov %*% t(jacobian))
6: FUN(X[[i]], ...)
5: lapply(seq_len(nrow(data)), function(datarow) {
       FUN <- gradient_factory(data = data[datarow, ], model = model, 
           variables = variables, type = type, weights = weights, 
           eps = eps, varslist = varslist, ...)
       jacobian <- jacobian(FUN, coef(model)[names(coef(model)) %in% 
           c("(Intercept)", colnames(vcov))], weights = weights, 
           eps = eps)
       diag(jacobian %*% vcov %*% t(jacobian))
   })
4: do.call("rbind", lapply(seq_len(nrow(data)), function(datarow) {
       FUN <- gradient_factory(data = data[datarow, ], model = model, 
           variables = variables, type = type, weights = weights, 
           eps = eps, varslist = varslist, ...)
       jacobian <- jacobian(FUN, coef(model)[names(coef(model)) %in% 
           c("(Intercept)", colnames(vcov))], weights = weights, 
           eps = eps)
       diag(jacobian %*% vcov %*% t(jacobian))
   }))
3: build_margins(model = model, data = data_list[[i]], variables = variables, 
       type = NULL, vcov = vcov, vce = vce, iterations = iterations, 
       unit_ses = unit_ses, eps = eps, varslist = varslist, ...)
2: margins.polr(mod, unit_ses = TRUE)
1: margins(mod, unit_ses = TRUE)

[vincentarelbundock]

For reference, the problem seems to be related to this line:

https://github.com/leeper/margins/blob/master/R/build_margins.R#L53

polr objects have several intercepts, so hard-coding "(Intercept)" won't work. In addition, polr objects store intercept estimates in object$zeta, whereas build_margins tries to feed coef(model) to jacobian. Note that the polr intercepts do appear in vcov(model).