i have 2 logistic , 2 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 y2 ordinal 1 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 marginal effects (mem approach) logistic models, used mfx package:
require(mfx) <- logitmfx(model, data=data, atmean=true) b <- logitmfx(modelinteraction, data=data, atmean=true) to calculate marginal effects ordered logistic models, used erer package:
require(erer) c <- ocme(mod) d <- ocme(modinteraction) what want is:
- plot results (i.e. variables)
a, b, c, , d. - show result 1 variable:
x1c(0,1) -- vary x1 between 0 , 1 -- while others hold @ mean ( both logistic , ordered logistic models).
the plot or table want create should this: figure 1 
the y axis in figure 1 indicates parameter estimate , x axis indicates name of variables
- i want plot interaction term in
b,d(i.e.x1*x5) figure similar this: figure 2
the y axis in figure 2 indicates probability difference , x axis minimum , maximum value of x5 (i.e. -10 +10)
i have been looking around solutions haven't been able find any. appreciate suggestions!
a reproducible sample (originally http://www.ats.ucla.edu/stat/data/binary.csv; i've made changes in order make more similar 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, create additional variable interaction term, because ocme() not appear behave when interaction term specified in formula.
data$x1_5 <- data$x1 * data$x5
then, fit model a, b, c, , d above (having changed x1*x5 in formulas 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 <- logitmfx(model, data=data, atmean=true) b <- logitmfx(modelinteraction, data=data, atmean=true) data$y2 <- as.factor(data$y2) # make y2 ordinal 1 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 can plot. make dataframe, out, contains coordinates want plot (the marginal effects , confidence intervals), based on logitmfx , ocme outputs. use 1.96 approximation critical levels, may or may not appropriate depending on size of 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 model b or est <- a$mfxest["x1",,drop=false] if want see estimates 1 variable. process similar ordered models, because marginal effects estimated each level of outcome variable, need plot level-specific marginal effects. estimated effects found in $out element of model fit, can put plotting code above loop, small modifications:
par(mfrow=c(1,3)) lvl <- 0 (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 bit more complicated, particularly confidence intervals. intuitive way (in mind) estimate uncertainty intervals using bootstrap (see king, tomz, , wittenberg 2000 in ajps, p. 352). uncertainty comes resampling data replacement. can write function bootstrapping, resample data , re-fit model:
bootstrap <- function(data, model) { newdata <- data[sample(rownames(data), nrow(data), replace=true),] fit <- polr(formula(model), data=newdata, method="logistic") } we fit model many times, each time using newly resampled dataset:
sims <- 1000 coefs <- replicate(sims, bootstrap(data, mod)) now have 1000 sets of parameter estimates. we'll use predict function generate new probabilities outcome variable. set 2 dataframes, x2, x3, , x4 take on mean values in data, x5 ranges -10 10 in 0.1 increments, , x1 0 , 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 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 can calculate marginal effects subtracting probabilities when x1=0 x1=1:
diffs <- lapply(1:sims, function(s) out_1[[s]] - out_0[[s]]) calculate means , 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, can plot 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, expected given estimates insignificant. interaction model, modify data_x1_0 , data_x1_1 take account interaction term (i.e. make new variable along lines of data_x1_0$x1_5 <- data_x1_0$x1 * data_x1_0$x5 -- zeros, , likewise data_x1_1), , modify coefs <- replicate(sims, bootstrap(data, mod)) use modinteraction instead of mod.
Comments
Post a Comment