1

我正在使用margins包（vignette）来很好地计算关于序数变量的边距。margins 包试图“移植 Stata（封闭源）边距的功能”。该视频（大约 4:25）显示，对于 Stata 中的序数概率模型，在我的示例中，我可以评估变量在序数变量的不同值下的边际效应x2。

我试过polr_1st_margins <- summary(margins(fit.polr, at = list(x2= 2:4)))了，这适用于示例数据。出于某种奇怪的原因，当我在我的实际数据上运行它时，我得到了一个错误。在这两种情况下x2都是一个因子变量（否则polr不会运行）。

我得到错误：Error in dat[, not_numeric, drop = FALSE] : incorrect number of dimensions即使变量在那里。

有人知道会发生什么吗？

示例代码：

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

编辑

我决定添加我的实际代码示例（也发布在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

不知何故，它没有找到变量，即使我可以看到它在那里。

数据

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

1 回答 1

0

我相信模型的margins代码中存在错误。polr你可能想试试这个marginaleffects包。 （免责声明：我是维护者。）

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
于 2021-11-22T01:50:06.597 回答