Как построить предельные эффекты (MEM) в R?

У меня есть две модели логистической регрессии и две модели упорядоченной логистической регрессии:

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) - измените X1 между 0 и 1 - в то время как другие сохранят свое среднее значение (как для логистической, так и для упорядоченной логистической модели).

График или таблица, которую я хочу создать, должны выглядеть так: Рисунок 1

Ось y на рисунке 1 указывает оценку параметра, а ось x указывает имя переменных.

Я также хочу изобразить термин взаимодействия в b и d (т.е. X1*X5), чтобы получить цифру, подобную этой: Рисунок 2

Ось y на рисунке 2 указывает разность вероятностей, а ось 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")

FKG 21.07.2016 источник

comment

Итак, чтобы использовать mod<- polr(X3 ~X1+X2+X4+X1*X4 ,data=data, Hess = TRUE), кажется, что X3 должен быть упорядоченным фактором (его нет в структуре, которую вы предоставили). Это намеренно, потому что в противном случае я предполагаю, что трансформирую его в единое целое. - Drey 25.07.2016

comment

Не могли бы вы уточнить, что именно вы подразумеваете под моим сюжетом все результаты в 1.? Потому что на данный момент кажется, что для a и b вы можете использовать plot(a$fit), plot(b$fit). Или вы хотите каким-то образом показать оценки от $mfxest (что и что?). Кроме того, что вы хотите показать от d? А что показывают оси в ваших примерах. (И исправьте код, который, к сожалению, неверен) - Drey 25.07.2016

comment

Привет @Drey Его надо преобразовать в порядковый, правильно. Под всеми результатами я имел в виду все переменные. У меня сейчас нет доступа к компьютеру, но через 2-3 часа он будет. Я смогу предоставить вам дополнительную информацию. Спасибо за ваше время! - FKG 25.07.2016

comment

@Drey Еще раз спасибо, что нашли время. Да, я хочу построить оценки из $mfxest и ocME для всех переменных в моделях в моделях a, b, c и d. Ось y на рисунке 1 указывает оценку параметра, а ось x указывает имена переменных. Ось y на рисунке 2 указывает разность вероятностей, а ось x указывает минимальное и максимальное значение переменной X4. Я понял, что образец не самый лучший, сразу поставлю образец получше. - FKG 25.07.2016

Ответы (2)

arrow_upward
3
arrow_downward

Я продолжу по частям.

Во-первых, я создаю дополнительную переменную для термина взаимодействия, потому что 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))

Назначьте est <- b$mfxest для модели b или 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 немного сложнее, особенно доверительные интервалы. На мой взгляд, наиболее интуитивно понятный способ оценки интервалов неопределенности - это использование начальной загрузки (см. King, Tomz, и Wittenberg 2000 в AJPS, стр. 352). Неопределенность возникает из-за повторной выборки данных с заменой. Мы можем написать функцию для выполнения начальной загрузки, в которой мы передискретизируем данные, а затем повторно подгоним модель:

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=0 из X1=1:

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.

Weihuang Wong 29.07.2016

comment

Большое спасибо, дорогой Вэйхуанг Вонг, за всю проделанную работу и за очень полезные R-коды. Чуть позже я смогу протестировать все это на моем реальном наборе данных. У меня могут возникнуть дополнительные вопросы - надеюсь, все в порядке. В любом случае, еще раз спасибо, что нашли время. Сразу отмечу твой ответ. - FKG; 30.07.2016

comment

Еще раз спасибо - работает очень хорошо. Для уточнения: (1) График взаимодействия не работает для логистической модели. Это из-за этой строки в функции начальной загрузки: fit <- **polr**(formula(model), data=newdata, method="logistic") Я получил ошибку, что ответ DV должен быть фактором. (2) `est‹ - a $ mfxest [X1`` drop = FALSE] `возможно ли отобразить более одной переменной? Я поигрался, но безрезультатно .. - FKG; 01.08.2016

comment

Также: используя эту строку для моих реальных данных data_means <- colMeans(data[,grep("X", names(data))], na.rm=TRUE) I get this error: object of type 'closure' is not subsettable. какова здесь функция X? Конечно, мои настоящие переменные имеют разные имена. Спасибо! - FKG; 01.08.2016

comment

(1) Да, вам следует изменить эту строку, чтобы использовать glm(...) для ваших биномиальных моделей (как в ваших моделях a и b). (2) Вы должны уметь делать что-то вроде est <- a$mfxest[c("X1","X2"),]. (3) grep("X", names(data)) возвращает индексы (позиции) имен столбцов, в которых есть X (например, X1, X2 и т. Д.). Таким образом, он исключает переменные Y. Но оглядываясь назад на то, как был написан код, простое выполнение colMeans(data, na.rm=TRUE) даст вам те же результаты, поэтому grep(...) не имеет значения. - Weihuang Wong; 01.08.2016

arrow_upward
1
arrow_downward

Если вы укажете свой термин взаимодействия таким образом, то ocME будет работать нормально:

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

Nikolos Gurney 19.09.2016

Как построить предельные эффекты (MEM) в R?

Ответы (2)

Похожие вопросы