Search code examples
rplotstatisticsglmmarginal-effects

How to plot marginal effects (MEM) in R?

I have two logistic and two ordered logistic regression models:

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

To calculate the marginal effects (MEM approach) for logistic models, I used the mfx package:

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

To calculate the marginal effects for ordered logistic models, I used the erer package:

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

What I want to do now is:

  1. plot all the results (i.e. all variables) for a, b, c, and d.
  2. show the result just for one variable: X1c(0,1) -- vary X1 between 0 and 1 -- while others hold at their mean ( for both logistic and ordered logistic models).

The plot or the table that I want to create should look like this: Figure 1 /Users/mac/Desktop/Skärmavbild 2016-07-21 kl. 22.59.47.png

The y axis in Figure 1 indicates parameter estimate and the x axis indicates the name of the variables

  1. I also want to plot the interaction term in b and d (i.e. X1*X5) to get a figure similar to this: Figure 2

enter image description here

The y axis in Figure 2 indicates probability difference and the x axis the minimum and maximum value of the X5 (i.e. -10 to +10)

I have been looking around for the solutions but haven't been able to find any. I would really appreciate any suggestions!

A reproducible sample (originally from http://www.ats.ucla.edu/stat/data/binary.csv; I've made some changes in order to make it more similar to my dataset): 

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

  • I'll proceed in parts.

    First, I create an additional variable for the interaction term, because ocME() does not appear to behave well when the interaction term is specified in the formula.

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

    Then, fit model a, b, c, and d as above (having changed X1*X5 in the formulas to 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)

    Now we can plot. I make a dataframe, out, that contains the coordinates that I want to plot (the marginal effects and the confidence intervals), based on the logitmfx and ocME outputs. I use 1.96 as an approximation for the critical levels, which may or may not be appropriate depending on the size of your dataset.

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

    Marginal effects 1

    Assign est <- b$mfxest for model b or est <- a$mfxest["X1",,drop=FALSE] if you just want to see the estimates for one variable. The process is similar for the ordered models, but because marginal effects are estimated for each level of the outcome variable, we need to plot level-specific marginal effects. The estimated effects are found in the $out element of the model fit, so we can put the plotting code above into a loop, with small modifications:

    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
}

    Marginal effects 2

    Figure 2 is a bit more complicated, particularly the confidence intervals. The most intuitive way (in my mind) to estimate the uncertainty intervals is using the bootstrap (see King, Tomz, and Wittenberg 2000 in AJPS, p. 352). The uncertainty comes from resampling the data with replacement. We can write a function to do the bootstrapping, where we resample the data and then re-fit the model:

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

    We fit the model many times, each time using a newly resampled dataset:

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

    Now we have 1000 sets of parameter estimates. We'll use the predict function to generate new probabilities for the outcome variable. We set up two dataframes, where X2, X3, and X4 take on the mean values in the data, X5 ranges from -10 to 10 in 0.1 increments, and X1 is 0 and 1 respectively.

    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

    Then use predict to get predicted probabilities:

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

    Now we can calculate the marginal effects by subtracting probabilities when X1=0 from X1=1:

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

    Calculate the means and the 95% interval:

    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)

    Finally, we can plot the results:

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

    Marginal effects 3

    Very unexciting, but to be expected given that the estimates are insignificant. To do this for the interaction model, modify data_X1_0 and data_X1_1 to take into account the interaction term (i.e. make a new variable along the lines of data_X1_0$X1_5 <- data_X1_0$X1 * data_X1_0$X5 -- which will be all zeros, and likewise for data_X1_1), and modify coefs <- replicate(sims, bootstrap(data, mod)) to use modInteraction instead of mod.