Как построить прогнозируемые поля, если они указаны с помощью «at»?

Question

Мы можем получить предельные эффекты линейной модели с помощью margins::margins() и можем выбрать интересующие переменные с помощью опции variables.

fit <- lm(mpg ~ factor(vs) + gear:factor(vs) + qsec, mtcars)

library(margins)
marg1 <- margins(fit, variables="vs")

> summary(marg1)
 factor    AME     SE      z      p   lower   upper
    vs1 4.8023 2.6769 1.7940 0.0728 -0.4443 10.0490

В пакете реализован метод plot.margins, поэтому мы можемГрафик предельных эффектов

plot(marg1)

at позволяет нам указать значения, по которым рассчитываются предельные эффекты:

marg2 <- margins(fit, variables="vs", at=list(gear=c(3, 4, 5)))

> summary(marg2)
 factor   gear    AME     SE      z      p   lower   upper
    vs1 3.0000 2.8606 3.3642 0.8503 0.3952 -3.7332  9.4544
    vs1 4.0000 5.6849 2.6713 2.1282 0.0333  0.4493 10.9206
    vs1 5.0000 8.5093 3.8523 2.2089 0.0272  0.9588 16.0597

Однако попытка построить эти указанные поля приведет к ошибке:

plot(marg2)

Error in `[.data.frame`(summ, , names(attributes(x)[["at"]]), drop = FALSE) : 
  undefined columns selected

Поскольку пакет margins претендует на то, что он является "R-портом полей Stata"'command' , я ожидаю, что график, похожий на тот, который дает Stata:

Итак, как мы можем построить прогнозируемые полякогда они указаны с помощью at?

edit:

Обратите внимание, что это не совсем обычный график взаимодействия, поскольку

with(mtcars[mtcars$gear %in% c(3, 4, 5), ], 
     interaction.plot(gear, vs, mpg, pch=rep(1, 2), type="b"))

дает другой вывод:

Answer 1

Ошибка возникает из-за ошибки в методе plot для объектов класса "margins", plot.margins.
Это попытка ее исправить.Изменения находятся в теле функции, просто выполните это или сохраните в файле "plotmargins.R", а затем source("plotmargins.R").

plot.margins <-
function (x, pos = seq_along(marginal_effects(x, with_at = FALSE)), 
    which = colnames(marginal_effects(x, with_at = FALSE)), labels = gsub("^dydx_", 
        "", which), horizontal = FALSE, xlab = "", ylab = "Average Marginal Effect", 
    level = 0.95, pch = 21, points.col = "black", points.bg = "black", 
    las = 1, cex = 1, lwd = 2, zeroline = TRUE, zero.col = "gray", 
    ...) 
{
    pars <- list(...)
    summ <- summary(x, level = level, by_factor = TRUE)
    MEs <- summ[, "AME", drop = TRUE]
    lower <- summ[, ncol(summ) - 1L]
    upper <- summ[, ncol(summ)]
    r <- max(upper) - min(lower)

    #--- changes start here
    nms <- intersect(names(summ), names(attributes(x)[["at"]]))
    at_levels <- unique(summ[, nms, drop = FALSE])
    #--- changes end here

    n_at_levels <- nrow(at_levels)
    if (n_at_levels > 1) {
        pos2 <- rep(pos, each = n_at_levels)
        pos2 <- pos2 + seq(from = -0.2, to = 0.2, length.out = n_at_levels)
    }
    else {
        pos2 <- pos
    }
    if (isTRUE(horizontal)) {
        xlim <- if ("xlim" %in% names(pars)) 
            xlim
        else c(min(lower) - 0.04 * r, max(upper) + 0.04 * r)
        ylim <- if ("ylim" %in% names(pars)) 
            xlim
        else c(min(pos2) - (0.04 * min(pos2)), max(pos2) + (0.04 * 
            max(pos2)))
    }
    else {
        xlim <- if ("xlim" %in% names(pars)) 
            xlim
        else c(min(pos2) - (0.04 * min(pos2)), max(pos2) + (0.04 * 
            max(pos2)))
        ylim <- if ("ylim" %in% names(pars)) 
            xlim
        else c(min(lower) - 0.04 * r, max(upper) + 0.04 * r)
    }
    if (isTRUE(horizontal)) {
        plot(NA, xlim = xlim, ylim = ylim, yaxt = "n", xlab = ylab, 
            ylab = xlab, las = las, ...)
        if (isTRUE(zeroline)) {
            abline(v = 0, col = zero.col)
        }
        points(MEs, pos2, col = points.col, bg = points.bg, pch = pch)
        axis(2, at = pos, labels = as.character(labels), las = las)
        mapply(function(pos, upper, lower, lwd) {
            segments(upper, pos, lower, pos, col = points.col, 
                lwd = lwd)
        }, pos2, upper, lower, seq(max(lwd), 0.25, length.out = length(MEs)))
    }
    else {
        plot(NA, xlim = xlim, ylim = ylim, xaxt = "n", xlab = xlab, 
            ylab = ylab, las = las, ...)
        if (isTRUE(zeroline)) {
            abline(h = 0, col = zero.col)
        }
        points(pos2, MEs, col = points.col, bg = points.bg, pch = pch)
        axis(1, at = pos, labels = as.character(labels), las = las)
        mapply(function(pos, upper, lower, lwd) {
            segments(pos, upper, pos, lower, col = points.col, 
                lwd = lwd)
        }, pos2, upper, lower, seq(max(lwd), 0.25, length.out = length(MEs)))
    }
    invisible(x)
}

Теперь ваш код и график.

source("plotmargins.R")

marg2 <- margins(fit, variables = "vs", 
                 at = list(gear = c(3, 4, 5)))

plot(marg2)