如何在 R 中绘制边际效应 (MEM)？

Question

我有两个逻辑回归模型和两个有序逻辑回归模型：

model <- glm(Y1 ~ X1+X2+X3+X4+X5, data = data, family = "binomial") #logistic
modelInteraction <- glm(Y1 ~ X1+X2+X3+X4+X5+X1*X5, data = data, family = "binomial") #logistic



        require(MASS)
        data$Y2 <- as.factor(data$Y2) # make the Y2 into a ordinal one 

mod<- polr(Y2 ~X1+X2+X3+X4+X5 ,data=data, Hess = TRUE) #ordered logistic
modInteraction<- polr(Y2~X1+X2+X3+X4+X5+X1*X5 ,data=data, Hess = TRUE) #ordered logistic

为了计算逻辑模型的边际效应（MEM 方法），我使用了 mfx 包：

require(mfx)
a <- logitmfx(model, data=data, atmean=TRUE)
b <- logitmfx(modelInteraction, data=data, atmean=TRUE)

为了计算有序逻辑模型的边际效应，我使用了 erer 包：

require(erer)    
c <- ocME(mod)
d <- ocME(modInteraction)

我现在要做的是：

绘制 a, b, c, and d 的所有结果（即所有变量）。
仅显示一个变量的结果：X1c(0,1)——在 0 和 1 之间改变 X1——而其他变量保持其均值（对于逻辑模型和有序逻辑模型）。

我想要创建的情节或 table 应该如下所示：图一

图1中的y轴表示参数估计，x轴表示变量名称

我还想绘制 b 和 d 中的交互项（即 X1*X5）以获得与此类似的图形：图 2

图2中的y轴表示概率差异，x轴表示X5的最小值和最大值（即-10到+10）

我一直在四处寻找解决方案，但一直找不到。如果有任何建议，我将不胜感激！

一个可重现的样本（最初来自http://www.ats.ucla.edu/stat/data/binary.csv；我做了一些更改以使其更类似于我的数据集）：

  > dput(data)
structure(list(Y1 = c(0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 
1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, NA, 1L, 0L, 1L, 
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 
0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 
1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 
0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 
0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, NA, 
0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 
0L, NA, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 
1L, 0L, 0L, 0L, 0L, 0L), Y2 = structure(c(1L, 3L, 2L, 2L, 1L, 
2L, 2L, 1L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 
2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 
1L, NA, 3L, 1L, 2L, 2L, 1L, 1L, 3L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 3L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 
1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 3L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 2L, 1L, 3L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 
3L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 3L, 1L, 1L, 1L, 
1L, 3L, 3L, 3L, 1L, 1L, 3L, 2L, 1L, 2L, 1L, 3L, 1L, 1L, 2L, 1L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, NA, 1L, 1L, 1L, 3L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 3L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 2L, 1L, 1L, 1L, 3L, 1L, 
3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 3L, 1L, 2L, 2L, 1L, 
1L, 3L, 1L, 2L, 2L, 1L, NA, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 2L, 
2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 
1L, 3L, 2L, 3L, 2L, 3L, 1L, 1L, 1L, 1L, 1L), .Label = c("0", 
"1", "2"), class = "factor"), X1 = c(0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 
1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), X2 = c(380L, 660L, 800L, 
640L, 520L, 760L, 560L, 400L, 540L, 700L, 800L, 440L, 760L, 700L, 
700L, 480L, 780L, 360L, 800L, 540L, 500L, 660L, 600L, 680L, 760L, 
800L, 620L, 520L, 780L, 520L, 540L, 760L, 600L, 800L, 360L, 400L, 
580L, 520L, NA, 520L, 560L, 580L, 600L, 500L, 700L, 460L, 580L, 
500L, 440L, 400L, 640L, 440L, 740L, 680L, 660L, 740L, 560L, 380L, 
400L, 600L, 620L, 560L, 640L, 680L, 580L, 600L, 740L, 620L, 580L, 
800L, 640L, 300L, 480L, 580L, 720L, 720L, 560L, 800L, 540L, 620L, 
700L, 620L, 500L, 380L, 500L, 520L, 600L, 600L, 700L, 660L, 700L, 
720L, 800L, 580L, 660L, 660L, 640L, 480L, 700L, 400L, 340L, 580L, 
380L, 540L, 660L, 740L, 700L, 480L, 400L, 480L, 680L, 420L, 360L, 
600L, 720L, 620L, 440L, 700L, 800L, 340L, 520L, 480L, 520L, 500L, 
720L, 540L, 600L, 740L, 540L, 460L, 620L, 640L, 580L, 500L, 560L, 
500L, 560L, 700L, 620L, 600L, 640L, 700L, 620L, 580L, 580L, 380L, 
480L, 560L, 480L, 740L, 800L, 400L, 640L, 580L, 620L, 580L, 560L, 
480L, 660L, 700L, 600L, 640L, 700L, 520L, 580L, 700L, 440L, 720L, 
500L, 600L, 400L, 540L, 680L, 800L, 500L, 620L, 520L, 620L, 620L, 
300L, 620L, 500L, 700L, 540L, 500L, 800L, 560L, 580L, 560L, 500L, 
640L, 800L, 640L, 380L, 600L, 560L, 660L, 400L, 600L, 580L, 800L, 
580L, 700L, 420L, 600L, 780L, 740L, 640L, 540L, 580L, 740L, 580L, 
460L, 640L, 600L, 660L, 340L, 460L, 460L, 560L, 540L, 680L, 480L, 
800L, 800L, 720L, 620L, 540L, 480L, 720L, 580L, 600L, 380L, 420L, 
800L, 620L, 660L, 480L, 500L, 700L, 440L, 520L, 680L, 620L, 540L, 
800L, 680L, 440L, 680L, 640L, 660L, 620L, 520L, 540L, 740L, 640L, 
520L, 620L, 520L, 640L, 680L, 440L, 520L, 620L, 520L, 380L, 560L, 
600L, 680L, 500L, 640L, 540L, 680L, 660L, 520L, 600L, 460L, 580L, 
680L, 660L, 660L, 360L, 660L, 520L, 440L, 600L, 800L, 660L, 800L, 
420L, 620L, 800L, 680L, 800L, 480L, 520L, 560L, NA, 540L, 720L, 
640L, 660L, 400L, 680L, 220L, 580L, 540L, 580L, 540L, 440L, 560L, 
660L, 660L, 520L, 540L, 300L, 340L, 780L, 480L, 540L, 460L, 460L, 
500L, 420L, 520L, 680L, 680L, 560L, 580L, 500L, 740L, 660L, 420L, 
560L, 460L, 620L, 520L, 620L, 540L, 660L, 500L, 560L, 500L, 580L, 
520L, 500L, 600L, 580L, 400L, 620L, 780L, 620L, 580L, 700L, 540L, 
760L, 700L, 720L, 560L, 720L, 520L, 540L, 680L, NA, 560L, 480L, 
460L, 620L, 580L, 800L, 540L, 680L, 680L, 620L, 560L, 560L, 620L, 
800L, 640L, 540L, 700L, 540L, 540L, 660L, 480L, 420L, 740L, 580L, 
640L, 640L, 800L, 660L, 600L, 620L, 460L, 620L, 560L, 460L, 700L, 
600L), X3 = c(3.61, 3.67, 4, 3.19, 2.93, 3, 2.98, 3.08, 3.39, 
3.92, 4, 3.22, 4, 3.08, 4, 3.44, 3.87, 2.56, 3.75, 3.81, 3.17, 
3.63, 2.82, 3.19, 3.35, 3.66, 3.61, 3.74, 3.22, 3.29, 3.78, 3.35, 
3.4, 4, 3.14, 3.05, 3.25, 2.9, NA, 2.68, 2.42, 3.32, 3.15, 3.31, 
2.94, 3.45, 3.46, 2.97, 2.48, 3.35, 3.86, 3.13, 3.37, 3.27, 3.34, 
4, 3.19, 2.94, 3.65, 2.82, 3.18, 3.32, 3.67, 3.85, 4, 3.59, 3.62, 
3.3, 3.69, 3.73, 4, 2.92, 3.39, 4, 3.45, 4, 3.36, 4, 3.12, 4, 
2.9, 3.07, 2.71, 2.91, 3.6, 2.98, 3.32, 3.48, 3.28, 4, 3.83, 
3.64, 3.9, 2.93, 3.44, 3.33, 3.52, 3.57, 2.88, 3.31, 3.15, 3.57, 
3.33, 3.94, 3.95, 2.97, 3.56, 3.13, 2.93, 3.45, 3.08, 3.41, 3, 
3.22, 3.84, 3.99, 3.45, 3.72, 3.7, 2.92, 3.74, 2.67, 2.85, 2.98, 
3.88, 3.38, 3.54, 3.74, 3.19, 3.15, 3.17, 2.79, 3.4, 3.08, 2.95, 
3.57, 3.33, 4, 3.4, 3.58, 3.93, 3.52, 3.94, 3.4, 3.4, 3.43, 3.4, 
2.71, 2.91, 3.31, 3.74, 3.38, 3.94, 3.46, 3.69, 2.86, 2.52, 3.58, 
3.49, 3.82, 3.13, 3.5, 3.56, 2.73, 3.3, 4, 3.24, 3.77, 4, 3.62, 
3.51, 2.81, 3.48, 3.43, 3.53, 3.37, 2.62, 3.23, 3.33, 3.01, 3.78, 
3.88, 4, 3.84, 2.79, 3.6, 3.61, 2.88, 3.07, 3.35, 2.94, 3.54, 
3.76, 3.59, 3.47, 3.59, 3.07, 3.23, 3.63, 3.77, 3.31, 3.2, 4, 
3.92, 3.89, 3.8, 3.54, 3.63, 3.16, 3.5, 3.34, 3.02, 2.87, 3.38, 
3.56, 2.91, 2.9, 3.64, 2.98, 3.59, 3.28, 3.99, 3.02, 3.47, 2.9, 
3.5, 3.58, 3.02, 3.43, 3.42, 3.29, 3.28, 3.38, 2.67, 3.53, 3.05, 
3.49, 4, 2.86, 3.45, 2.76, 3.81, 2.96, 3.22, 3.04, 3.91, 3.34, 
3.17, 3.64, 3.73, 3.31, 3.21, 4, 3.55, 3.52, 3.35, 3.3, 3.95, 
3.51, 3.81, 3.11, 3.15, 3.19, 3.95, 3.9, 3.34, 3.24, 3.64, 3.46, 
2.81, 3.95, 3.33, 3.67, 3.32, 3.12, 2.98, 3.77, 3.58, 3, 3.14, 
3.94, 3.27, 3.45, 3.1, 3.39, 3.31, 3.22, 3.7, 3.15, 2.26, 3.45, 
2.78, 3.7, 3.97, 2.55, 3.25, 3.16, NA, 3.5, 3.4, 3.3, 3.6, 3.15, 
3.98, 2.83, 3.46, 3.17, 3.51, 3.13, 2.98, 4, 3.67, 3.77, 3.65, 
3.46, 2.84, 3, 3.63, 3.71, 3.28, 3.14, 3.58, 3.01, 2.69, 2.7, 
3.9, 3.31, 3.48, 3.34, 2.93, 4, 3.59, 2.96, 3.43, 3.64, 3.71, 
3.15, 3.09, 3.2, 3.47, 3.23, 2.65, 3.95, 3.06, 3.35, 3.03, 3.35, 
3.8, 3.36, 2.85, 4, 3.43, 3.12, 3.52, 3.78, 2.81, 3.27, 3.31, 
3.69, 3.94, 4, 3.49, 3.14, NA, 3.36, 2.78, 2.93, 3.63, 4, 3.89, 
3.77, 3.76, 2.42, 3.37, 3.78, 3.49, 3.63, 4, 3.12, 2.7, 3.65, 
3.49, 3.51, 4, 2.62, 3.02, 3.86, 3.36, 3.17, 3.51, 3.05, 3.88, 
3.38, 3.75, 3.99, 4, 3.04, 2.63, 3.65, 3.89), X4 = c(3L, 3L, 
1L, 4L, 4L, 2L, 1L, 2L, 3L, 2L, 4L, 1L, 1L, 2L, 1L, 3L, 4L, 3L, 
2L, 1L, 3L, 2L, 4L, 4L, 2L, 1L, 1L, 4L, 2L, 1L, 4L, 3L, 3L, 3L, 
1L, 2L, 1L, 3L, NA, 3L, 2L, 2L, 2L, 3L, 2L, 3L, 2L, 4L, 4L, 3L, 
3L, 4L, 4L, 2L, 3L, 3L, 3L, 3L, 2L, 4L, 2L, 4L, 3L, 3L, 3L, 2L, 
4L, 1L, 1L, 1L, 3L, 4L, 4L, 2L, 4L, 3L, 3L, 3L, 1L, 1L, 4L, 2L, 
2L, 4L, 3L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 4L, 2L, 
2L, 3L, 3L, 3L, 4L, 3L, 2L, 2L, 1L, 2L, 3L, 2L, 4L, 4L, 3L, 1L, 
3L, 3L, 2L, 2L, 1L, 3L, 2L, 2L, 3L, 3L, 3L, 4L, 1L, 4L, 2L, 4L, 
2L, 2L, 2L, 3L, 2L, 3L, 4L, 3L, 2L, 1L, 2L, 4L, 4L, 3L, 4L, 3L, 
2L, 3L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 2L, 1L, 2L, 3L, 2L, 2L, 
2L, 2L, 2L, 1L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 4L, 2L, 2L, 3L, 
3L, 3L, 3L, 4L, 2L, 2L, 4L, 2L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 4L, 
2L, 2L, 3L, 4L, 3L, 4L, 3L, 2L, 1L, 4L, 1L, 3L, 1L, 1L, 3L, 2L, 
4L, 2L, 2L, 3L, 2L, 3L, 1L, 1L, 1L, 2L, 3L, 3L, 1L, 3L, 2L, 3L, 
2L, 4L, 2L, 2L, 4L, 3L, 2L, 3L, 1L, 2L, 2L, 2L, 4L, 3L, 2L, 1L, 
3L, 2L, 1L, 3L, 2L, 2L, 3L, 3L, 4L, 4L, 2L, 4L, 4L, 3L, 2L, 3L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 3L, 2L, 3L, 2L, 3L, 2L, 1L, 
2L, 2L, 3L, 1L, 4L, 2L, 2L, 3L, 4L, 4L, 2L, 4L, 1L, 4L, 4L, 4L, 
2L, 2L, 2L, 1L, 1L, 3L, 1L, NA, 2L, 3L, 2L, 3L, 2L, 2L, 3L, 4L, 
1L, 2L, 2L, 3L, 3L, 2L, 3L, 4L, 4L, 2L, 2L, 4L, 4L, 1L, 3L, 2L, 
4L, 2L, 3L, 1L, 2L, 2L, 2L, 4L, 3L, 3L, 1L, 3L, 3L, 1L, 3L, 4L, 
1L, 3L, 4L, 3L, 4L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 2L, 
2L, 1L, 2L, 1L, 3L, 3L, 1L, 1L, 2L, NA, 1L, 3L, 3L, 3L, 1L, 2L, 
2L, 3L, 1L, 1L, 2L, 4L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 3L, 2L, 3L, 2L, 3L, 2L, 2L, 3L), X5 = c(10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, 
-7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, -7L, 
-7L, -7L, -6L, 7L, -7L, -7L, -7L, 7L, 7L, 7L, 7L, 7L, 2L, -2L, 
-2L, -2L, -2L, 0L, 3L, 5L, 5L, 5L, 5L, 0L, 0L, 6L, 6L, 6L, 6L, 
6L, 5L, 5L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 10L, 10L, 10L, 10L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 0L, 0L, 0L, 0L, 0L, 4L, 4L, 4L, 6L, 
6L, 6L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 6L, 6L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
8L, 8L, 8L, -1L, 6L, 6L, 6L, 6L, 6L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, -3L, -3L, -3L, -3L, -3L, 
-3L, -3L, -3L, -3L, -3L, -3L, -3L, -3L, -3L, -4L, -4L, -4L, -4L, 
-4L, -4L, -4L, -4L, -4L, -4L, -4L, -4L, -4L, -4L, -4L, -5L, -5L, 
-5L, -5L, -5L, -5L, -5L, -5L, -5L, -5L, -5L, -5L, -5L, -8L, -8L, 
-8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, -8L, 
-9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, 
-9L, -9L, -10L, -10L, -10L, -10L, -10L, -10L, -10L, -10L, -10L, 
-10L, -10L, -10L, -10L)), .Names = c("Y1", "Y2", "X1", "X2", 
"X3", "X4", "X5"), row.names = c(NA, -400L), class = "data.frame")

Answer 1

我将分部分进行。

首先，我为交互项创建了一个附加变量，因为当在公式中指定交互项时，ocME() 似乎表现不佳。

data$X1_5 <- data$X1 * data$X5

然后，拟合模型 a、b、c 和 d（已将公式中的 X1*X5 更改为 X1_5).

require(MASS)  # polr
require(mfx)   # logitmfx
require(erer)  # ocME

model <- glm(Y1 ~ X1+X2+X3+X4+X5, data = data, family = "binomial") #logistic
modelInteraction <- glm(Y1 ~ X1+X2+X3+X4+X5+X1_5, data = data, family = "binomial") #logistic
a <- logitmfx(model, data=data, atmean=TRUE)
b <- logitmfx(modelInteraction, data=data, atmean=TRUE)

data$Y2 <- as.factor(data$Y2) # make the Y2 into a ordinal one 
mod<- polr(Y2 ~X1+X2+X3+X4+X5 ,data=data, Hess = TRUE) #ordered logistic
modInteraction<- polr(Y2~X1+X2+X3+X4+X5+X1_5 ,data=data, Hess = TRUE) #ordered logistic
c <- ocME(mod)
d <- ocME(modInteraction)

现在我们可以绘图了。我制作了一个数据框 out，其中包含我想要绘制的坐标（边际效应和置信区间），基于 logitmfx 和 ocME 输出。我使用 1.96 作为临界水平的近似值，这可能合适也可能不合适，具体取决于您的数据集的大小。

est <- a$mfxest
par(mfrow=c(1,1))
out <- data.frame(mean=est[,1],
                  lower=est[,1]-1.96*est[,2],
                  upper=est[,1]+1.96*est[,2])
plot(x=1:nrow(out), y=out$mean, ylim=c(min(out$lower), max(out$upper)), 
    xaxt="n", ylab="Marginal effects", xlab="", las=2)
abline(h=0, col="grey")
arrows(x0=1:nrow(out), y0=out$lower, y1=out$upper, code=3, angle=90, length=.05)
axis(1, at=1:nrow(out), labels=rownames(out))

为模型b分配est <- b$mfxest，如果您只想查看一个变量的估计值，则分配est <- a$mfxest["X1",,drop=FALSE]。有序模型的过程类似，但由于边际效应是针对结果变量的每个水平估算的，我们需要绘制特定水平的边际效应。估计效果在模型拟合的 $out 元素中找到，因此我们可以将上面的绘图代码放入循环中，稍作修改：

par(mfrow=c(1,3))
lvl <- 0
for (est in c$out[1:3]) {
    out <- data.frame(mean=est[,1],
                      lower=est[,1]-1.96*est[,2],
                      upper=est[,1]+1.96*est[,2])
    plot(x=1:nrow(out), y=out$mean, 
        ylim=c(min(out$lower), max(out$upper)),
        xlim=c(.5, nrow(out)+.5),
        xaxt="n", ylab="", xlab="", las=2,
        main=paste("Marginal effects on Level", lvl))
    abline(h=0, col="grey")
    arrows(x0=1:nrow(out), y0=out$lower, y1=out$upper, code=3, angle=90, length=.05)
    axis(1, at=1:nrow(out), labels=rownames(out))
    lvl <- lvl + 1
}

图 2 有点复杂，尤其是置信区间。估计不确定性区间的最直观方法（在我看来）是使用 bootstrap（参见 AJPS 中的 King, Tomz, and Wittenberg 2000，第 352 页）。不确定性来自于通过替换对数据进行重新采样。我们可以编写一个函数来执行 bootstrapping，我们在其中重新采样数据，然后重新拟合模型：

bootstrap <- function(data, model) {
    newdata <- data[sample(rownames(data), nrow(data), replace=TRUE),]
    fit <- polr(formula(model), data=newdata, method="logistic")
}

我们多次拟合模型，每次都使用新的重采样数据集：

sims <- 1000
coefs <- replicate(sims, bootstrap(data, mod))

现在我们有1000组参数估计。我们将使用 predict 函数为结果变量生成新的概率。我们设置了两个数据框，其中 X2、X3 和 X4 取数据中的平均值，X5 范围从 -10 到 10，增量为 0.1，并且X1分别为0和1。

data_means <- colMeans(data[,grep("X", names(data))], na.rm=TRUE)
data_X1_0 <- data.frame(X1=0,
                        X2=data_means["X2"],
                        X3=data_means["X3"],
                        X4=data_means["X4"],
                        X5=seq(-10, 10, .1))
data_X1_1 <- data_X1_0
data_X1_1$X1 <- 1

然后用predict得到预测概率：

out_0 <- lapply(coefs, function(fit) predict(fit, data_X1_0, type="probs"))
out_1 <- lapply(coefs, function(fit) predict(fit, data_X1_1, type="probs"))

现在我们可以通过从X1=1中减去X1=0时的概率来计算边际效应：

diffs <- lapply(1:sims, function(s) out_1[[s]] - out_0[[s]])

计算均值和 95% 区间：

diffs <- array(unlist(diffs), 
    dim = c(nrow(diffs[[1]]), ncol(diffs[[1]]), length(diffs)))
means <- apply(diffs, MARGIN=c(1,2), mean)
upper <- apply(diffs, MARGIN=c(1,2), quantile, .975)
lower <- apply(diffs, MARGIN=c(1,2), quantile, .025)

最后，我们可以绘制结果：

for (i in 1:3) {
    plot(x=seq(-10, 10, .1), y=means[,i], type="l", 
        ylim=c(min(lower[,i]), max(upper[,i])), xlab="", ylab="")
    lines(x=seq(-10, 10, .1), y=upper[,i], lty=2)
    lines(x=seq(-10, 10, .1), y=lower[,i], lty=2)   
}

非常平淡，但考虑到估计值微不足道，这是可以预料的。要为交互模型执行此操作，请修改 data_X1_0 和 data_X1_1 以考虑交互项（即，沿 data_X1_0$X1_5 <- data_X1_0$X1 * data_X1_0$X5 行创建一个新变量——这将全为零， data_X1_1) 也是如此，并修改 coefs <- replicate(sims, bootstrap(data, mod)) 以使用 modInteraction 而不是 mod.

Answer 2

如果你这样指定你的交互项，那么 ocME 会很好地播放：

# Ordered logistic
modInteraction <- polr(Y2~X1+X2+X3+X4+X5+I(X1*X5), data=data, Hess=TRUE)

如何在 R 中绘制边际效应 (MEM)？

How to plot marginal effects (MEM) in R?

r

plot

statistics

glm

marginal-effects