FKG - 1 year ago 209
R Question

# 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:
`X1`
c(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

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

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

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

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

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

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`.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download