How to use margins package to evaluate marginal affects at different values of the dependent variable

Question

I am using the margins package (vignette) to well, calculate margins, with respect to an ordinal variable. The margins package is an attempt to "port the functionality of Stata’s (closed source) margins". This video (around 4:25) shows that for an ordinal probit model in Stata, I can evaluate the marginal effect of a variable at different values of the ordinal variable, in my example x2.

I tried polr_1st_margins <- summary(margins(fit.polr, at = list(x2= 2:4))) and this works on the example data. For some weird reason, when I run it on my actual data I get an error. In both cases x2 is a factor variable (otherwise polr won't run).

I get the error: Error in dat[, not_numeric, drop = FALSE] : incorrect number of dimensions even when the variable is there.

Does someone have an idea what could be going on?

Example code:

library(sure) # for residual function and sample data sets
library(MASS) # for polr function
library(margins)
rm(df1)
df1 <- df1
df1$x2 <- df2$y
df1$x3 <- df2$x
df1$y <- df3$x/10
fit.polr <- polr(x2 ~ x + x3, data = df1, method = "probit")
summary(margins(fit.polr))

polr_1st_margins <- summary(margins(fit.polr, at = list(x2= 2:4)))

EDIT

I have decided to add a sample of my actual code (also posted on GitHub).

fit.polr2 <- polr(ordinal_dep_var ~ Dummy + Continuous + Dummy2 + as.factor(industry) + Urbanisation_Dummy + Size_Dummy, data = df2, method = "probit", Hess=TRUE)
summary(margins(fit.polr2))

polr_1st_margins <- summary(margins(fit.polr2, at = list(ordinal_dep_var= 0:3)))

Error in dat[, not_numeric, drop = FALSE] : 
  incorrect number of dimensions

Somehow, it does not find the variable, even though I can see that it is there.

DATA

df2 <- structure(list(ordinal_dep_var = structure(c(4L, 3L, 4L, 1L, 
3L, 1L, 1L, 3L, 1L, 4L, 4L, 4L, 1L, 3L, 4L, 2L, 4L, 2L, 2L, 1L, 
2L, 1L, 1L, 1L, 3L, 4L, 3L, 2L, 2L, 1L, 1L, 2L, 4L, 2L, 1L, 4L, 
3L, 1L, 2L, 3L, 4L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 2L, 4L, 2L, 1L, 
4L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 3L, 
2L, 3L, 3L, 1L, 4L, 4L, 4L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 2L, 4L, 3L, 3L, 1L, 1L, 1L, 3L, 2L, 3L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 4L, 3L, 2L, 2L, 3L, 1L, 1L, 1L, 
1L, 3L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("0", 
"1", "2", "3"), class = c("ordered", "factor")), Dummy = c(0, 
0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0), Continuous = c(15.6862745098039, 
41.6666666666667, 26.9230769230769, 14.8514851485149, 32.1428571428571, 
20, 43.75, 0, 0, 0, 80, 100, 21.4285714285714, 23.8095238095238, 
28.125, 0, 30.3030303030303, 25, 100, 100, 100, 13.3333333333333, 
66.6666666666667, 33.3333333333333, 55.5555555555556, 72.2222222222222, 
57.3033707865169, 17.7777777777778, 47.6190476190476, 17.7777777777778, 
41.6666666666667, 40, 20, 8.33333333333333, 16.6666666666667, 
40, 100, 0, 50, 0, 0, 7.69230769230769, 0, 0, 0, 0, 0, 0, 0, 
33.3333333333333, 20, 1.84089414858646, 1.84089414858646, 1.84089414858646, 
30, 20, 0, 33.3333333333333, 33.3333333333333, 0, 0, 0, 100, 
50, 22.2222222222222, 0, 0, 50, 50, 46.1538461538462, 44.4444444444444, 
0, 5.55555555555556, 0, 0, 47.3684210526316, 43.75, 18.1818181818182, 
0, 42.8571428571429, 14.2857142857143, 0, 50, 0, 0, 50, 20, 50, 
100, 0, 42.8571428571429, 20, 25, 33.3333333333333, 0, 0, 0, 
66.6666666666667, 0, 25.9259259259259, 0, 33.3333333333333, 0, 
100, 25, 0, 0, 10, 100, 50, 33.3333333333333, 75, 0, 0, 40, 0, 
33.3333333333333, 28.5714285714286, 0, 0, 28.5714285714286, 0, 
28.5714285714286, 0, 0, 9.09090909090909, 30, 66.6666666666667, 
50, 0, 20, 50, 0, 33.3333333333333, 0, 66.6666666666667, 18.1818181818182, 
28.5714285714286, 36.9230769230769, 14.2857142857143, 36.3636363636364, 
0, 0, 7.69230769230769), Dummy2 = structure(c(0, 0, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 
0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 
1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 
0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 
0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 
1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0), label = "Gift/informal payment requested: tax inspectorate?", format.stata = "%14.0f", class = c("haven_labelled", 
"vctrs_vctr", "double"), labels = c(Yes = 1, No = 2)), industry = structure(c(3, 
7, 2, 7, 19, 17, 9, 7, 9, 11, 12, 10, 5, 5, 3, 12, 3, 5, 4, 1, 
1, 11, 5, 7, 9, 9, 7, 9, 5, 9, 4, 7, 9, 11, 11, 22, 12, 21, 19, 
11, 10, 7, 19, 23, 20, 20, 21, 24, 19, 1, 21, 21, 21, 21, 21, 
21, 6, 6, 21, 3, 17, 19, 20, 12, 21, 20, 12, 10, 10, 21, 21, 
19, 3, 21, 21, 10, 10, 21, 5, 20, 21, 19, 22, 15, 6, 21, 21, 
10, 19, 21, 20, 3, 6, 10, 24, 24, 21, 23, 21, 21, 11, 21, 3, 
3, 21, 24, 21, 1, 10, 22, 19, 17, 3, 20, 23, 1, 22, 21, 21, 23, 
21, 23, 21, 22, 22, 21, 10, 13, 10, 15, 10, 24, 24, 7, 10, 10, 
22, 21, 21, 23, 21, 22, 7, 21), label = "Industry", format.stata = "%34.0g", class = c("haven_labelled", 
"vctrs_vctr", "double"), labels = c(Textiles = 1, Leather = 2, 
Garments = 3, Agroindustry = 4, Food = 5, Beverages = 6, `Metals and machinery` = 7, 
Electronics = 8, `Chemicals and pharmaceutics` = 9, Construction = 10, 
`Wood and furniture` = 11, `Non-metallic and plastic materials` = 12, 
Paper = 13, `Sport goods` = 14, `IT services` = 15, `Other manufacturing` = 16, 
Telecommunications = 17, `Accounting and finance` = 18, `Advertising and marketing` = 19, 
`Other services` = 20, `Retail and wholesale trade` = 21, `Hotels and restaurants` = 22, 
Transport = 23, `Real estate and rental services` = 24, `Mining and quarrying` = 25, 
`Auto and auto components` = 26, `Other transport equipment` = 27, 
`Other unclassified` = 99)), Urbanisation_Dummy = structure(c(2L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 
3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 
1L, 3L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 3L, 1L, 1L, 2L, 3L, 
3L, 3L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 2L, 1L, 1L, 1L, 2L, 3L, 2L, 
3L, 3L, 3L, 3L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 1L, 2L, 3L, 2L, 2L, 3L, 1L, 2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 
3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 3L, 2L, 1L, 1L, 2L, 
2L, 1L, 1L, 3L, 1L, 3L, 2L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 1L, 3L, 1L, 1L, 3L, 3L, 1L, 3L, 3L, 2L, 3L, 1L, 1L), .Label = c("City > 250", 
"50k-250k", "< 50k"), class = "factor"), Size_Dummy = structure(c(3L, 
3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 3L, 3L, 
3L, 1L, 3L, 3L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 
3L, 3L, 3L, 3L, 3L, 2L, 3L, 1L, 2L, 3L, 3L, 2L, 1L, 1L, 1L, 3L, 
3L, 2L, 1L, 1L, 1L, 1L, 1L, 3L, 2L, 3L, 3L, 3L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 3L, 3L, 1L, 3L, 3L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 
1L, 3L, 1L, 3L, 2L, 1L, 1L, 1L, 1L, 1L, 3L, 2L, 3L, 3L, 1L, 3L, 
1L, 1L, 1L, 3L, 3L, 2L, 3L, 1L, 2L, 1L, 3L, 1L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 1L, 3L, 3L, 3L, 
2L, 3L, 3L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L), .Label = c("Employees: < 10", 
"Employees: 10-19", "Employees: 20+"), class = "factor")), row.names = c(NA, 
-144L), class = c("data.table", "data.frame"))

Answer 1

I believe that there is a bug in the margins code for polr models. You might want to try the marginaleffects package instead. (Disclaimer: I am the maintainer.)

library(MASS)
library(marginaleffects)

mod <- polr(factor(gear) ~ hp, data = mtcars, method = "probit")

mfx <- marginaleffects(mod, type = "probs")

summary(mfx)

## Average marginal effects 
##    type Group Term     Effect Std. Error z value Pr(>|z|)     2.5 %    97.5 %
## 1 probs     3   hp  0.0009916  0.0011329  0.8753  0.38142 -0.001229 0.0032122
## 2 probs     4   hp -0.0003951  0.0004637 -0.8521  0.39415 -0.001304 0.0005137
## 3 probs     5   hp -0.0005965  0.0007168 -0.8322  0.40530 -0.002001 0.0008084
## 
## Model type:  polr 
## Prediction type:  probs